多充液可機動柔性航天器模塊化建模和耦合動力學分析

摘 要：針對現(xiàn)代航天器因攜帶多類型充液貯箱及可機動柔性附件而導致的動力學建模困難，文中基于參數(shù)辨識理論和有限元方法實現(xiàn)了該類復(fù)雜航天器的模塊化建模。首先，考慮航天器運行時貯箱內(nèi)液體的小幅晃動問題，利用球形貯箱小幅晃動的參數(shù)化模型推導出球形充液貯箱小幅晃動的動力學方程；其次，根據(jù)薄板小幅振動理論，運用有限單元法建立柔性附件四邊形板單元動力學模型，將附件相對航天器變化的姿態(tài)角代入坐標轉(zhuǎn)換矩陣，構(gòu)造柔性附件作大范圍機動時的時變坐標轉(zhuǎn)換矩陣。再次，基于凱恩方程推導航天器整體系統(tǒng)的動力學狀態(tài)方程，利用MATLAB軟件編制出相應(yīng)的模塊化建模程序。最后，通過數(shù)值算例分析，研究典型構(gòu)型航天器中柔性附件以不同方式機動的系統(tǒng)整體耦合動力學性能，驗證了該建模方法的通用性、適用性和準確性。

關(guān)鍵詞：充液柔性航天器；參數(shù)化模型；有限單元法；凱恩方程；模塊化建模

中圖分類號：V412.4；V448.2 DOI：10.16375/j.cnki.cn45-1395/t.2024.01.001

0 引言

現(xiàn)代航天器不斷向多任務(wù)、模塊化和長壽命方向發(fā)展，攜帶諸多設(shè)備和大尺寸柔性附件的航天器是復(fù)雜的剛-液-柔耦合系統(tǒng)，準確建立系統(tǒng)動力學模型不僅能描述航天器系統(tǒng)的耦合動力學特性，也能為系統(tǒng)控制律的設(shè)計提供理論依據(jù)和前期研究基礎(chǔ)。

航天器的典型構(gòu)型通常為中心剛體攜帶多個充液貯箱，外側(cè)連接太陽能帆板。太陽能電池板因具有較大長寬尺度和極小厚度的幾何尺寸特征，根據(jù)具體長寬比可簡化為Euler-Bernoulli梁模型[1-3]或Kirchhoff-Love薄板模型[4]。附件相對主剛體運動的耦合效應(yīng)[3，6-7]通過定義相對運動廣義坐標[3]或構(gòu)造時變坐標轉(zhuǎn)換矩陣[6]等途徑對航天器動力學模型進行改良。現(xiàn)代航天器所攜帶的液體燃料占據(jù)總質(zhì)量的較大比重[8]，等效力學模型被廣泛應(yīng)用于液體的晃動效應(yīng)的簡化和替代，例如彈簧-質(zhì)量-阻尼模型[9-11]、單擺模型[12-13]和復(fù)合擺模型[14-15]。區(qū)別于傳統(tǒng)等效模型參數(shù)與機械結(jié)構(gòu)參數(shù)相對應(yīng)的限制，液體晃動參數(shù)化模型[16-18]為實驗或仿真數(shù)據(jù)逆向辨識所得到，適用性更強，而且對仿真數(shù)據(jù)的辨識精度可達92%以上[19]，具有較好的研究價值和應(yīng)用前景。對于非完整約束系統(tǒng)動力學方程的建立，Kane等[20-21]在推導系統(tǒng)動力學方程時，提出一種基于偏速度和偏角速度為系統(tǒng)自由度的建模方法，又稱為凱恩方程。相較于傳統(tǒng)拉格朗日法，凱恩方程步驟程序化、物理意義清楚，廣泛應(yīng)用于航天器系統(tǒng)動力學中[22]。

現(xiàn)階段，柔性附件充液航天器的動力學研究存在以下局限性：1）等效力學模型適用范圍存在一定局限，而參數(shù)化模型可以實時辨識并修改參數(shù)模型，描述更準確。2）拉格朗日方程推導多自由度耦合系統(tǒng)動力學方程較為困難且繁瑣。因此，本文基于文獻[19]的結(jié)論，構(gòu)造球形充液貯箱小幅振動參數(shù)化模型；運用有限單元法推導柔性附件動力學方程，根據(jù)附件相對航天器的時變姿態(tài)角計算時變坐標轉(zhuǎn)換矩陣；最終由凱恩方程模塊化地建立航天器系統(tǒng)的耦合動力學方程，研究系統(tǒng)狀態(tài)變量在柔性部件做不同機動和不同數(shù)量下的動力學響應(yīng)。

1 航天器系統(tǒng)剛-液-柔耦合動力學

方程推導

航天器系統(tǒng)的典型構(gòu)型通常簡化為中心剛體帶多個充液貯箱，并外接柔性附件，如圖1所示。定義以下參考系：慣性參考系N，[n1、n2、n3]為慣性參考系基矢量；與航天器本體固連的本體參考系B，其原點位于主剛體質(zhì)心O，[b1、b2、b3]為本體參考系的基矢量，初始時刻本體系基矢量與慣性系對應(yīng)基矢量方向一致；附件局部坐標系S，其原點位于附件在航天器的安裝點[OS]，[s11、s12、s13]為附件局部參考系的基矢量；附件有限元模型中第i個四邊形板單元的局部坐標系[Ki]，其原點位于四邊形板單元幾何中心，[ki1][、][ki2][、][ki3]為單元參考系基矢量，單元參考系基矢量與附件局部參考系S對應(yīng)基矢量方向一致。

1.1 航天器主體姿態(tài)動力學建模

航天器的軌道機動由質(zhì)心O相對于慣性系N的平動進行描述，姿態(tài)調(diào)整由本體系B相對于慣性系N的轉(zhuǎn)動進行描述。定義[rO=NrTON=rO1，rO2，rO3T]

[n1，n2，n3]為主剛體質(zhì)心O相對慣性系的矢徑，其中列矩陣[NrO=rO1，rO2，rO3]為矢徑[rO]投影到慣性坐標系N在其基矢量方向上的分量陣，左上角標N表示矢量向慣性參考系投影；[N=n1，n2，n3]為慣性系N的基矢量組合。同理，矢徑[rO]可由本體系B基矢量表示為：[rO=BrTOB]，其中[BrTO、B]的定義與上述N同理。航天器在軌運行時其相對慣性系的矢徑快速變化，在研究小擾動對航天器的影響時，將本體參考系下計算矢徑響應(yīng)結(jié)果即為相對穩(wěn)定運行的影響。對主剛體質(zhì)心O相對慣性系N的絕對速度[vO]由矢徑[rO]對時間求導得到，并投影到本體坐標系B中，則有

[vO=drOdt=BvTOB+Bω×0BrOTB，] （1）

式中：[Bω×0]為角速度在本體系投影的分量矩陣[Bω0]的反對稱矩陣，其中[ω0]為本體參考系B相對慣性參考系N旋轉(zhuǎn)的角速度。主剛體的絕對加速度[aO]可通過對速度求導后得到，即

[aO=BaTOB+Bω×0BvOTB+Bα×0BrOTB+]

[Bω×0Bω×0BrOTB，] （2）

式中：[Bα×0]為角速度在本體系投影的分量陣[Bα0]的反對稱矩陣，其中[α0]為本體參考系B相對慣性參考系N旋轉(zhuǎn)的角加速度。

1.2 柔性附件動力學建模

考慮三軸姿態(tài)穩(wěn)定航天器在軌運行時擾動較小，帆板變形基于小變形理論描述，則無需考慮材料大幅變形時的動力剛化效應(yīng)。構(gòu)造太陽能帆板四邊形薄板有限單元模型，將太陽能帆板劃分為多個較小四邊形板單元并對節(jié)點編號。根據(jù)有限元法，板單元各節(jié)點具有3個自由度，即單元變形時節(jié)點產(chǎn)生3個方向的位移：相對[ki1、ki2]的角位移[θ1、θ2]以及垂直于板平面的撓度w。撓度僅與該點在單元局部參考系的平面坐標相關(guān)，即[w（ki1，ki2）]；根據(jù)Kirchhoff-Love板理論，薄板厚度方向各點具有相同撓度變形w。

考慮第[sl]個柔性附件上一點P，該點在慣性參考系的矢徑可由矢量疊加法得到，即

[rP=rO+rsl+rslP+uP=BrTOB+BrTslB+]

[CBSSrTslP+SuTPB，] （3）

式中：[rsl=BrTslB]為第[sl]個柔性附件安裝點[OSl]相對主剛體質(zhì)心的矢徑，并在本體參考系B下表示，[Brsl]為[rsl]在本體參考系B投影的分量矩陣；[rslP=CBSSrTslPB]為P點相對附件安裝點的矢徑，[SrslP]為[rslP]在附件局部參考系S投影的分量矩陣，[CBS]為附件局部參考系與本體系的坐標轉(zhuǎn)換矩陣，轉(zhuǎn)換關(guān)系為[S=CBSB]；[uP=CBSSuTPB]為P點發(fā)生變形時相對初始位置的位移矢量，列矩陣[SuP]為[uP]在附件局部參考系S投影的分量矩陣。根據(jù)有限元法，第i個四邊形板單元上一點P垂直于板平面的撓度變形w由第i個單元的節(jié)點位移[δei]和該點在單元局部參考系坐標[ki1，ki2]共同決定，即[SuP=0，0，wki1， ki2T=0， 0， Nki1， ki2?δeiT]，其中[Nx， y]為4節(jié)點板單元形函數(shù)。同理，質(zhì)點P在慣性參考系的矢徑[rP]對時間求導得到絕對速度[vP]，并將其投影到本體坐標系后，有

[vP=drPdt=BvToB+Bω×0BroTB+Bω×0BrslTB+]

[CBSSω×0+Sω×1SrP+SuPTB+CBSSuTPB，] （4）

式中：[Sω×0和Sω×1]分別為主剛體角速度和附件相對角速度在局部參考系S投影的分量矩陣[Sω0和Sω1]的反對稱矩陣，[ω1]為柔性附件相對主剛體的旋轉(zhuǎn)角速度。

柔性附件變形后質(zhì)點相對慣性系的絕對加速度[aP]可描述為

[aP=dvPdt=BaToB+Bω×0BvoTB+Bα×0BrOTB+][Bω×0Bω×0BrOTB+Bα×0BrslTB+Bω×0Bω×0BrslTB+][Sα×0+Sα×1SrslP+SuPTS+2Sω×0+Sω×1SuPTS+][Sω×0+Sω×1Sω×0+Sω×1SrslP+SuPTS+SuTPS，] （5）

式中：[α1]為柔性附件相對主剛體的旋轉(zhuǎn)角速度；[Sα×0和Sα×1]分別為主剛體角加速度和附件相對角加速度在局部參考系S投影的分量矩陣[Sα0和Sα1]的反對稱矩陣。

有限單元法中將作用于板單元的分布載荷由虛功原理等效為單元節(jié)點載荷，由式（5）質(zhì)點加速度[aP]和形函數(shù)[Nx，y]在單元內(nèi)積分得到單元節(jié)點載荷矩陣：[Fe=NTρhaPdS]。按節(jié)點編號將單元節(jié)點載荷矩陣組裝得到附件總體載荷矩陣[Fsl]，總體質(zhì)量矩陣M和總體剛度矩陣K[23]。定義[δsl]為附件節(jié)點位移列矩陣，聯(lián)立得到柔性附件四邊形板單元振動變形的動力學方程為

[Mδsl+Cδsl+Kδsl=Fsl.] （6）

1.3 液體晃動等效力學模型

航天器液體燃料貯箱安裝位置如圖1所示。由于三軸姿態(tài)穩(wěn)定航天器運行時產(chǎn)生的擾動較小，航天器控制系統(tǒng)消耗的燃料忽略不計，即假設(shè)液體燃料質(zhì)量不變。由文獻[19]的結(jié)論，基于勢流理論和傅里葉-貝塞爾函數(shù)展開系數(shù)條件，在忽略高階項時球形貯箱液體晃動狀態(tài)方程如式（7）所示。顯然該狀態(tài)方程的結(jié)構(gòu)與傳統(tǒng)物理等效力學模型動力學方程相似，例如單擺等效力學模型等。

[ddtb1（t）a1（t）=01-ω2-2ωζb1（t）a1（t）+0Aax（t），] （7）

式中：[b1（t）， a1（t）T]為系統(tǒng)狀態(tài)向量；[ζ]為一階模態(tài)的阻尼比；[ω]為系統(tǒng)固有頻率；[A]為輸入矩陣系數(shù)；[ax（t）]為橫向晃動加速度。

在研究液體晃動效應(yīng)時大多不需要推導流體域各處情況，通常將晃動液體視為整體參與晃動，因此根據(jù)狀態(tài)方程，設(shè)參與晃動的液體等效質(zhì)心m在外激勵下做小幅振動時，質(zhì)心相對貯箱球心以固定長度晃動；未參與晃動的靜止液體的等效質(zhì)心設(shè)為[m1]。由矢量疊加法，質(zhì)心m在慣性系中的矢徑如下

[rm=rO+rsP+rl+um=BrTOB+BrTsPB+BrTlB+BuTmB，]

（8）

式中：[rsP=BrTsPB]為球形貯箱球心相對主剛體質(zhì)心的矢徑，列矩陣[BrsP]為矢徑[rsP]在本體參考系B投影的分量矩陣；[rl=BrTlB]為充液貯箱靜置時質(zhì)心m相對球形貯箱球心的矢徑，列矩陣[Brl]為矢徑[rl]在本體參考系B投影的分量矩陣；[um=BuTmB]為質(zhì)心m晃動時偏離初始位置的位移矢量，列矩陣[Bum]是[um]在本體參考系B投影的分量矩陣。

在文獻[19]中，參數(shù)化模型對球形充液貯箱受橫向外激勵晃動的輸出響應(yīng)辨識效果較好，而且小幅晃動下質(zhì)心m擺動的角位移可近似等效為沿[b1、b2]方向的線位移，即[Bum=[um1，um2，0]T]。同理，晃動液體等效質(zhì)點m在慣性系下的速度和加速度通過對時間求導后，得到

[vm=drmdt=BvTOB+Bω×0BrsP+Brl+BumTB+BuTmB，] （9）

[am=dvmdt=BaTOB+Bω×0BvOTB+][Bα×0BrsP+Brl+BumTB+Bω×0Bω×0BrsP+Brl+BumTB+][Bω×0BumTB+BuTmB.] （10）

最終，球形充液貯箱小幅晃動（包括待辨識參數(shù)）的等效動力學狀態(tài)方程可整理為

[m00mum1um2+c00cum1um2=-mBam（1）-mBam（2）]， （11）

式中：c為液體晃動等效阻尼系數(shù)；[Bam（1）、Bam（2）]為質(zhì)心m的加速度本體參考系的投影在[b1、b2]方向上的分量。

2 基于凱恩方程的航天器耦合系統(tǒng)

動力學方程推導

凱恩將偏速度和偏角速度作為系統(tǒng)自由度，并聯(lián)立廣義主動力和廣義慣性力在偏速度和偏角速度方向的投影，最終得到系統(tǒng)動力學方程。

2.1 偏速度和偏角速度

偏速度和偏角速度統(tǒng)一用[v（r）i]表示，其中右上角標（r）表示第r個偏速度。對于航天器主剛體相對慣性系的平動和轉(zhuǎn)動，定義偏速度和偏角速度：

[vji=BvO（j）?bj， " " " " "j=1，2，3，] （12）

[v3+ji=Bω0（j）?bj]， [j=1，2，3]. （13）

設(shè)航天器共安裝[sl]個附件，每個柔性附件劃分n個節(jié)點，四邊形板單元每個節(jié)點含3個自由度。對于第s個柔性附件的第j個節(jié)點的振動變形，進一步定義偏速度和偏角速度：

[v6+3n（s-1）+3（j-1）+1i=θ1j?b1，]

[s=1，2，…，sl;j=1，2，…，n]， （14）

[v6+3n（s-1）+3（j-1）+2i=θ2j?b2，]

[s=1，2，…，sl;j=1，2，…，n]， （15）

[v6+3n（s-1）+3（j-1）+3i=wj?b3，]

[s=1，2，…，sl;j=1，2，…，n.] （16）

另設(shè)航天器總共攜帶q個部分充液貯箱。對于第k個充液貯箱的參與晃動液體質(zhì)心m相對主剛體的運動，可定義如下偏速度和偏角速度：

[v6+3nsl+2（k-1）+1i=um1?b1， k=1，2，…，q]， （17）

[v6+3nsl+2（k-1）+2i=um2?b2][， k=1，2，…，q]. （18）

2.2 廣義慣性力和廣義主動力

對于第j個偏速度方向，航天器耦合系統(tǒng)所受的廣義慣性力[Fj]的計算公式為

[Fj=-mBBaO?v（j）i-IB?NαB0+NωB0×IB?NωB0?v（j）i-]

[s=1slSNaP?v（j）idmP-k=1qmkBamk?v（j）i，] （19）

式中：[mB]為主剛體的質(zhì)量；[IB]為主剛體關(guān)于本體系的二階慣性張量；[mP]為附件上任意質(zhì)點P的質(zhì)量；[mk]為第k個貯箱內(nèi)液體晃動等效質(zhì)心的質(zhì)量；[amk]為第k個充液貯箱晃動液體等效質(zhì)點的晃動加速度。

航天器系統(tǒng)所受的廣義作用力包括液體晃動的反作用力，以及柔性附件振動變形的彈性作用力。對于第j個偏速度方向，航天器耦合系統(tǒng)所受的廣義作用力[Fj*]的計算公式為

[Fj*=FsP?v（j）i+Fsl?v（j）i]， （20）

式中：[FsP=k=1qmk（um1b1+um2b2）+c（um1b1+um2b2）]為液體晃動等效質(zhì)心的作用力；[Fsl]為柔性附件振動變形對支座的反作用力，計算公式為[Fsl=s=1slCSBKj，：sδsls]，其中[CSB]為附件局部參考系S到本體系B的坐標轉(zhuǎn)換矩陣，（j，：）表示第j行所有列的全部元素。

2.3 基于凱恩方程的系統(tǒng)動力學方程建模

將相同偏速度方向下的廣義主動力和廣義慣性力分量聯(lián)立得到系統(tǒng)動力學方程。在第j個偏速度方向下的凱恩方程可寫為

[Fj=F*j]. （21）

將式（12）—式（18）中求得的偏速度代入式（19）和式（20），可建立描述系統(tǒng)運動的動力學方程組。將航天器系統(tǒng)所有自由度的響應(yīng)按偏速度的順序組合成位移列矩陣u，最后將得到的動力學方程組整理成如下形式，

[Mu+Cu+Ku=F]， （22）

式中：M、C、K分別為系統(tǒng)質(zhì)量矩陣、阻尼矩陣和剛度矩陣。質(zhì)量矩陣M由主剛體、柔性附件和液體晃動參數(shù)模型的質(zhì)量矩陣對角排列構(gòu)成，且包含了廣義主動力中系數(shù)為廣義速率導數(shù)的項；同理，阻尼矩陣C和剛度矩陣K的排列構(gòu)成類似，且包含了廣義主動力中系數(shù)為廣義速率的項和系數(shù)為自由度位移的項。考慮MATLAB微分方程ODE求解器僅適用于一階方程，需對系統(tǒng)方程改寫為等效一階方程組。系統(tǒng)廣義位移和廣義速度記為系統(tǒng)的狀態(tài)向量[u，uT]，將上述動力學方程降階為

[IOOMddtuu=OI-K-Cuu+OF]. （23）

3 數(shù)值仿真

對上述建立的含柔性附件充液柔性航天器的系統(tǒng)動力學方程，編寫相應(yīng)MATLAB程序。在算例1中，航天器各部件參數(shù)和構(gòu)型參照文獻[24]，航天器所受外激勵如下所示

[M02（t）10， 0lt;tlt;10， " " " " " " " " " " " " " -10， 20lt;tlt;30， " " " " " " " " " " 0， t≤0， 10≤t≤20， t≥30. ] （24）

編程計算航天器系統(tǒng)在外激勵作用下的耦合動力學響應(yīng)如圖2—圖5所示。其中，航天器速度與角速度響應(yīng)如圖2與圖3所示，柔性附件自由邊中點的位移響應(yīng)如圖4所示，液體晃動等效質(zhì)心位移響應(yīng)如圖5所示。

編程求解結(jié)果與文獻[24]相比較，由于忽略了帆板小幅振動時的阻尼效應(yīng)，在外激勵作用結(jié)束后帆板仍進行小幅的振動變形，而航天器響應(yīng)與附件位移響應(yīng)曲線的趨勢和數(shù)值大小與參考文獻結(jié)果相接近，初步驗證了推導和程序的準確性。

太陽能帆板屬于大撓度附件，當帆板相對航天器本體運動時對其姿態(tài)角產(chǎn)生較大影響。在算例2、3中，考慮柔性附件進行機動時，計算航天器動力學響應(yīng)。算例2的航天器構(gòu)型如圖1所示，初始時刻各部件處于靜止狀態(tài)。航天器系統(tǒng)包含4個球形充液貯箱，貯箱參數(shù)參考文獻[25]，其球心相對本體坐標系原點的坐標如表1所示。

航天器本體安裝一個機動太陽能帆板，該柔性附件安裝點為[rsl=[0，0.15，0]T]，安裝角為[θ01=[0，0，π 2]T]。柔性附件的尺寸和材料參數(shù)如表2所示。

機動附件的驅(qū)動類型采用Bang-Bang驅(qū)動[26]，柔性附件繞局部參考系的[s12]軸相對航天器旋轉(zhuǎn)。其相應(yīng)的角加速度的變化如下所示：

[α（t）=0.021， 0lt;tlt;5， " " " " " "-0.021， 5lt;tlt;10， " " " 0， t≤0， t=5， t≥10. " ] （25）

編程計算柔性附件機動時航天器動力學響應(yīng)如圖6—圖7所示，充液貯箱晃動質(zhì)心位移響應(yīng)如圖8—圖9所示，柔性附件自由邊中點位移響應(yīng)如圖10所示。

附件繞坐標軸[s12]進行旋轉(zhuǎn)機動時，振動變形的反作用力矩為本體參考系的[b1]方向，因此結(jié)合反作用力對質(zhì)心的力矩，在附件機動時角速度響應(yīng)[ω01]較大。附件的機動形式和驅(qū)動形式相關(guān)，振動趨勢與文獻[26]中算例相接近。由于部分充液貯箱等效晃動質(zhì)心一般并不位于球心，因此主剛體運動產(chǎn)生的牽連載荷各不相同，產(chǎn)生復(fù)雜的晃動現(xiàn)象。

算例3：柔性附件繞局部參考系的[s11]軸相對航天器旋轉(zhuǎn)，驅(qū)動和其余參數(shù)與算例2對應(yīng)參數(shù)相同。編程計算柔性附件機動時航天器動力學響應(yīng)如圖11—圖12所示，充液貯箱晃動質(zhì)心位移響應(yīng)如圖13—圖14所示，柔性附件自由邊中點位移響應(yīng)如圖15所示。

附件繞坐標軸[s11]進行旋轉(zhuǎn)機動時，振動變形的反作用力矩為本體參考系的[b2]方向，相較于前文算例，附件機動時產(chǎn)生一定的角速度響應(yīng)[ω02]。由于附件相對[s11]的慣性矩更小，因此在相同驅(qū)動下各部件的動力學響應(yīng)也更小。

4 結(jié)論

本文將三軸姿態(tài)穩(wěn)定充液柔性航天器為研究對象，考慮航天器運行過程中的小幅擾動，基于球形貯箱液體晃動參數(shù)化模型推導充液貯箱晃動的動力學方程，運用有限元法建立板類柔性附件的四邊形單元模型，基于Kirchhoff-Love板理論描述單元位移。由凱恩方程建立航天器剛-液-柔耦合系統(tǒng)動力學模型。仿真對比分析充液柔性航天器在單個板類附件沿不同轉(zhuǎn)軸機動下的動力學響應(yīng)，結(jié)果顯示航天器姿態(tài)角和軌道狀態(tài)參數(shù)、液體晃動效應(yīng)、附件機動與柔性變形的相互耦合作用，驗證了該方法建立的耦合系統(tǒng)動力學模型的準確性。體現(xiàn)了該方法對不同附件驅(qū)動類型、不同附件數(shù)量航天器的模塊化建模特性。

參考文獻

[1] GHORBANI H，VATANKHAH R，F(xiàn)ARID M. General planar motion modeling and control of a smart rigid-flexible satellite considering large deflections[J]. Nonli-near Dynamics，2022，108：911-939.

[2] GHORBANI H，VATANKHAH R，F(xiàn)ARID M. Adaptive nonsingular fast terminal sliding mode controller design for a smart flexible satellite in general planar motion[J]. Aerospace Science and Technology，2021，119. Doi：10.1016/j.ast.2021.107100.

[3] MANCINI M，CAPELLO E. Adaptive sliding mode-based control system for flexible spacecraft[J]. Proceedings of the American Control Conference，2021（2021-May）：2968-2973.

[4] DA FONSECA I M，RADE D A，GOES L C S，et al. Attitude and vibration control of a satellite containing flexible solar arrays by using reaction wheels，and piezoelectric transducers as sensors and actuators[J]. Acta Astronautica，2017（139）：357-366.

[5] HE G Q，CAO D Q，CAO Y T，et al. Dynamic modeling and orbit maneuvering response analysis for a three-axis attitude stabilized large scale flexible spacecraft installed with hinged solar arrays[J]. Mechanical Systems and Signal Processing，2022，162：108083.

[6] WANG J，LI D X，JIANG J P. First order coupled dynamic model of flexible space structures with time-varying configurations[J]. Acta Astronautica，2017，132：117-123.

[7] ZHANG G W，WANG F，CHEN J，et al. Fixed-time sliding mode attitude control of a flexible spacecraft with rotating appendages connected by magnetic bearing[J]. Mathematical Biosciences and Engineering，2022，19（3）：2286-2309.

[8] BURK T A. Cassini orbit trim maneuvers at saturn-overview of attitude control flight operations[M]. AIAA Guidance，Navigation，and Control Conference，2011.

[9] ALLARD C，RAMOS M D，SCHAUB H，et al. Spacecraft dynamics integrating hinged solar panels and lumped-mass fuel slosh model[C]//AIAA/AAS Astrodynamics Specialist Conference 2011，2011.

[10] CHIBA M，MAGATA H. Coupled pitching dynamics of flexible space structures with on-board liquid sloshing[J]. Acta Astronautica，2021，181：151-166.

[11] COLAGROSSI A，LAVAGNA M. Integrated vibration suppression attitude control for flexible spacecrafts with internal liquid sloshing[J]. Multibody System Dynamics，2020，51（2）：123-157.

[12] SUN Y，ZHOU D，WANG J D，et al. Liquid sloshing in a cylindrical tank with multiple baffles under horizontal and pitching motions[J]. International Journal of Applied Mechanics，2020，12（7）：190-197.

[13] GASBARRI P，SABATINI M，PISCULLI A. Dynamic modelling and stability parametric analysis of a flexible spacecraft with fuel slosh[J]. Acta Astronautica，2016，127：141-159.

[14] LIU F，YUE B Z，TANG Y，et al. 3DOF-rigid-pendulum analogy for nonlinear liquid slosh in spherical propellant tanks[J]. Journal of Sound and Vibration，2019，460：114907.

[15] LU Y，YUE B Z，MA B L. Improved moving pulsating ball equivalent model for large-amplitude liquid slosh[J]. AIAA Journal，2022，60（8）：5004-5011.

[16] 周鳳霞.帶貯液罐類重型車輛操縱動力學與穩(wěn)定性研究[D].柳州：廣西科技大學，2020.

[17] 黃中烈，吳文軍，周鳳霞，等.罐車動力學建模與剛-液耦合特性研究[J].現(xiàn)代制造工程，2020（8）：20-26.

[18] 李超，吳文軍，周鳳霞，等.重型罐車內(nèi)液體晃動參數(shù)化數(shù)學模型辨識[J].廣西科技大學學報，2021，32（2）：65-71.

[19] 王佐，吳文軍，姜丙超，等.球形貯箱內(nèi)液體橫向晃動問題的系統(tǒng)辨識與參數(shù)化等效建模[J].廣西科技大學學報，2022，33（3）：1-7.

[20] KANE T R，LEVINSON D A. Formulation of equations of motion for complex spacecraft[J]. Journal of Guidance，Control，and Dynamics，1980，3（2）：99-112.

[21] KANE T R，RYAN R R，BANERJEE A K. Dynamics of a cantilever beam attached to a moving base[J]. Journal of Guidance，Control，and Dynamics，1987，10（2）：139-151.

[22] LIU F，YUE B Z，BANERJEE A K，et al. Large motion dynamics of in-orbit flexible spacecraft with large-amplitude propellant slosh[J]. Journal of Guidance，Control，and Dynamics，2020，43（3）：438-450.

[23] 韓清凱，孫偉，王伯平，等.機械結(jié)構(gòu)有限單元法基礎(chǔ)[M].北京：科學出版社，2013.

[24] 閆玉龍，申云峰，岳寶增.含板類柔性附件的充液航天器動力學研究[J].宇航學報，2020，41（5）：531-540.

[25] 苗楠，李俊峰，王天舒.橫向激勵下液體大幅晃動建模分析[J].宇航學報，2016，37（3）：268-274.

[26] 吳文軍.帶多充液貯箱及機動柔性附件航天器耦合動力學研究[D].北京：北京理工大學，2015.

Modular modeling and coupled dynamical study for liquid-filling spacecraft with multiple maneuverable flexible appendages

LUO Houlin， WU Wenjun*， WANG Zuo

（School of Mechanical and Automotive Engineering， Guangxi University of Science and Technology，

Liuzhou 545616， China）

Abstract： This paper focuses on the problem of dynamic modeling for modern spacecraft with multi-type liquid-filling tanks and maneuverable flexible appendages. Based on parameter identification theory and finite element method， the dynamic modeling for complex spacecraft is deduced. Firstly， in view of small disturbance during spacecraft operation， the dynamic equations of small-amplitude liquid sloshing in spherical tank are derived by using the parametric model for liquid sloshing. Secondly， based on the theory for small deformation of plate， the dynamic quadrilateral plate element model of flexible appendage is established using finite element method. The time-varying coordinate transformation matrix of flexible appendage relative to spacecraft is constructed by substituting the attitude angle of the appendage into the transformation matrix. Then， the dynamic state equations of the spacecraft system are derived based on Kane's equation， and the corresponding modular modeling program is compiled by MATLAB software. Finally， through numerical example analysis， the coupling dynamic performance of the typical configuration spacecraft with flexible appendages maneuvering in different ways is studied， and the universality， availability and accuracy of the modeling method are confirmed.

Keywords： liquid-filling flexible spacecraft; parametric model; finite element method; Kane's equation; modular modeling

（責任編輯：黎 婭、于艷霞）