2.3.1 Spacecraft hub

The in-orbit motion of the spacecraft body comprises two forms: translational and rotational. The generalized speed of the spacecraft body undergoing rigid-body translation is defined as u _i (i = 1, 2, 3). As shown in Eq. (5). the translational velocity of the spacecraft body in the inertial reference frame is ^N v ^{B *} , as expressed in Eq. (6).

(5) \begin{align} u_{i}={}^{N}{\boldsymbol{v}}{_{}^{B^{*}}}\cdot \boldsymbol{b}_{i}\left(i=1,2,3\right) \\[-8pt] \nonumber \end{align}

(6) \begin{align} {}^{N}{\boldsymbol{v}}{_{}^{B^{*}}}=u_{i}\cdot \boldsymbol{b}_{i}\left(i=1,2,3\right) \\[12pt] \nonumber \end{align}

The generalized speed of rotation is u _{3 + i} (i = 1, 2, 3) As shown in Eq. (7), the angular velocity of the spacecraft body’s attitude rotation relative to the inertial reference frame is ^N $\boldsymbol{\omega}$ ^B As shown in Eq. (8).

(7) \begin{align} u_{3+i}={}^{N}{\boldsymbol{\omega }}{_{}^{B}}\cdot \boldsymbol{b}_{i}\left(i=1,2,3\right) \\[-8pt] \nonumber \end{align}

(8) \begin{align} {}^{N}{\boldsymbol{\omega }}{_{}^{B}}=u_{3+i}\cdot \boldsymbol{b}_{i}\left(i=1,2,3\right) \\[2pt] \nonumber \end{align}

2.3.2 Liquid-filled tank

As shown in Figure 1, the position vector p ^{NP _k} of the equivalent center of mass of the liquid in the k _th tank within the inertial reference frame is:

(9) \begin{align}\boldsymbol{p}^{N{P_{k}}}=\boldsymbol{p}^{N{B^{*}}}+\boldsymbol{p}^{{B^{*}}{T_{k}}}+\boldsymbol{p}^{{T_{k}}{P_{k}}} \end{align}

in which p ^{NB *} represents the position vector of the spacecraft’s center of mass relative to the inertial frame; p ^{B * T _k} denotes the position vector of the k _th tank’s centroid relative to the spacecraft’s center of mass, determined by the tank’s mounting position within the spacecraft and thus a constant vector, defined as Eq. (10); p ^{T _k P _k} signifies the position vector of the k _th tank’s liquid sloshing-centered equivalent mass point P _k relative to the k _th tank’s centroid, defined as Eq. (11).

(10) \begin{align} \boldsymbol{p}^{{B^{*}}{T_{k}}}=X_{k}\boldsymbol{b}_{1}+Y_{k}\boldsymbol{b}_{2}+Z_{k}\boldsymbol{b}_{3}\left(k=1\sim 2\right) \\[-8pt] \nonumber \end{align}

(11) \begin{align} \boldsymbol{p}^{{T_{k}}{P_{k}}}=x_{k}\boldsymbol{b}_{1}+y_{k}\boldsymbol{b}_{2}+z_{k}\boldsymbol{b}_{3}\left(k=1\sim 2\right) \\[10pt] \nonumber \end{align}

Let the generalized speed of motion of the liquid-filled tank’s concentrated equivalent center of mass relative to the inertial reference frame be u _{3 + 3k + i} (i = 1, 2, 3):

(12) \begin{align}u_{3+3k+i}=\left({{}^{B}{{\rm d}}{} \over {\rm d}t}\boldsymbol{p}^{{T_{k}}{P_{k}}}\right)\cdot \boldsymbol{b}_{i}{\enspace}\left(i=1,2,3; \, k=1,2\right) \end{align}

The velocity ^N v ^Pk of the equivalent center of mass P _k relative to the inertial reference frame can be obtained by differentiating p ^{NP _k} with respect to time, i.e., Eq. (13).

(13) \begin{align}\begin{array}{c}{}^{N}{\boldsymbol{v}} {_{}^{P_{k}}}= \dfrac{{}^{N}{{\rm d}}}{ {\rm d}t}\boldsymbol{p}^{N{P_{k}}}= \dfrac{{}^{N}{{\rm d}}} {{\rm d}t}\left(\boldsymbol{p}^{N{B^{*}}}+\boldsymbol{p}^{{B^{*}}{T_{k}}}+\boldsymbol{p}^{{T_{k}}{P_{k}}}\right)={}^{N}{\boldsymbol{v}}{_{}^{B^{*}}}+{}^{N}{\boldsymbol{\omega }}{_{}^{B}}\times \left(\boldsymbol{p}^{{B^{*}}{T_{k}}}+\boldsymbol{p}^{{T_{k}}{P_{k}}}\right)+ \dfrac{{}^{B}{{\rm d}}}{{\rm d}t} \boldsymbol{p}^{{T_{k}}{P_{k}}}\\[8pt] =\boldsymbol{b}_{1}[u_{1}+u_{5}(Z_{k}+z_{k})-u_{6}(Y_{k}+y_{k})+u_{4+3k}] \quad \quad\quad \quad\quad \quad\quad \quad\quad \quad\quad \quad\quad\\[5pt] +\boldsymbol{b}_{2}[u_{2}-u_{4}(Z_{k}+z_{k})+u_{6}(X_{k}+x_{k})+u_{5+3k}] \quad\quad \quad\quad \quad\quad \quad\quad \quad\quad \quad \\[5pt] +\boldsymbol{b}_{3}[u_{3}+u_{4}(Y_{k}+y_{k})-u_{5}(X_{k}+x_{k})+u_{6+3k}] \quad\quad \quad\quad \quad\quad \quad\quad \quad\quad \quad \end{array} \end{align}

2.3.3 Flexible rod

As shown in Figure 5 the rotation of the rod around any point R _n in space is modeled as successive rotations about two mutually perpendicular axes. This rotation can be decomposed into two Euler rotations: The local coordinate frame R _n : { r _n1, r _n2, r _n3} rotates α _n rad around the $\boldsymbol{r}^{0}_{n2}$ axis from its initial position aligned with $R_{n}^{0}$ : { $\boldsymbol{r}^{0}_{n1}$ , $\boldsymbol{r}^{0}_{n2}$ , $\boldsymbol{r}^{0}_{n3}$ }, then moves to the intermediate coordinate frame $R_n^{*}$ : { $\boldsymbol{r}^{*}_{n1}$ , $\boldsymbol{r}^{*}_{n2}$ , $\boldsymbol{r}^{*}_{n3}$ }. It then rotates $\beta$ _n rad around the $\boldsymbol{r}^{*}_{n3}$ axis to reach its current position, where $R_{n}^{0}$ : { $\boldsymbol{r}^{0}_{n1}$ , $\boldsymbol{r}^{0}_{n2}$ , $\boldsymbol{r}^{0}_{n3}$ } is the mounting coordinate frame; R _n : { r _n1, r _n2, r _n3} is the local coordinate frame. Let the transformation matrix converting the flexible rod from $R_{n}^{0}$ to R _n be denoted as, then ${}^{{R_{n}^{0}}} \boldsymbol{C}^{Rn}$ is:

Figure 5.

The n _th flexible rod rotates relative to the rigid body.

(14) \begin{align}{}^{R_{n}^{0}} {\boldsymbol{C}} ^{R_{n}}=\left[\begin{array}{ccc} \cos \beta _{n} &\quad \sin \beta _{n} &\quad 0\\[3pt] -\sin \beta _{n} &\quad \cos \beta _{n} &\quad 0\\[3pt] 0 &\quad 0 &\quad 1 \end{array}\right]\left[\begin{array}{ccc} \cos a_{n} &\quad 0 &\quad -\sin a_{n}\\[3pt] 0 &\quad 1 &\quad 0\\[3pt] \sin a_{n} &\quad 0 &\quad \cos a_{n} \end{array}\right] \end{align}

Based on this, let the coordinate transformation matrix from the body-fixed coordinate frame B: { b ₁, ${\boldsymbol{b}_{2}}$ , ${\boldsymbol{b}_{3}}$ } to the mounting coordinate frame $R_{n}^{0}$ : { $\boldsymbol{r}^{0}_{n1}$ , $\boldsymbol{r}^{0}_{n2}$ , $\boldsymbol{r}^{0}_{n3}$ } be ${}^{B} \boldsymbol{C}^{R_{n}^{0}}$ . Combining with Eq. (14) yields the coordinate transformation matrix from the body-fixed coordinate frame B: { ${\boldsymbol{b}_{1}}$ , ${\boldsymbol{b}_{2}}$ , ${\boldsymbol{b}_{3}}$ } to the local coordinate frame R _n : { r _n1, r _n2, r _n3} as:

(15) \begin{align}{}^{B}{\boldsymbol{C}}{_{}^{R_{n}}}={}^{R_{n}^{0}}{\boldsymbol{C}^{Rn}} \,{}^{B}{\boldsymbol{C}}{_{}^{R_{n}^{0}}} \end{align}

Thus, the angular velocity of the flexible appendage relative to the inertial reference frame is given by:

(16) \begin{align}{}^{N}{\boldsymbol{\omega }}{_{}^{R}}={}^{N}{\boldsymbol{\omega }}{_{}^{B}}+\dot{\alpha }_{n}{\boldsymbol{r}_{n{2}}^{0}}+\dot{\beta }_{n}\boldsymbol{r}_{n3}=\left(\begin{array}{l} u_{4}\\[3pt] u_{5}\\[3pt] u_{6} \end{array}\right)+{\left({}^{B}{\boldsymbol{C}}{_{}^{R_{n}^{0}}}\right)^{\!\!\!{\rm T}}\left(\begin{array}{ccc} \cos a_{n} &\quad {\rm 0} &\quad -\sin a_{n}\\[3pt] {\rm 0} &\quad {\rm 1} &\quad {\rm 0}\\[3pt] \sin a_{n} &\quad {\rm 0} &\quad \cos a_{n} \end{array}\right)}^{\!\!\!\!\!{\rm T}}\left(\begin{array}{l} {\rm 0}\\[3pt] {\rm 0}\\[3pt] \dot{\beta } \end{array}\right) \end{align}

where $\dot{\alpha }_{n}={{}^{R_{n}^{0}}{{\rm d}}{} \over {\rm d}t}\alpha _{n}, \dot{\beta }_{n}={{}^{R}{{\rm d}}{} \over {\rm d}t}\beta _{n}$

In this paper, we employ the geometrically nonlinear theory based on the assumptions of Euler–Bernoulli beam theory to describe the structural vibrations of the flexible appendages. Assuming the large-scale flexible rod to be a uniform-section homogeneous cantilever beam, we incorporate the axial contraction induced by lateral deformation into the analysis. Consider an arbitrary mass element d $m_{n}^{\textrm{rod}}$ on the flexible rod. When the appendage is undeformed, its position point is H _n . After deformation, its position point becomes $H_{n}^{\textrm{def}}$ , with a deformation displacement of $\boldsymbol{d}_{n}^{\textrm{rod}}$ , as shown in Figure 6.

Figure 6.

Schematic diagram of flexible rod deformation.

Therefore, after elastic deformation of the flexible rod, the position of any mass element d $m_{n}^{\textrm{rod}}$ in the inertial reference frame can be expressed as:

(17) \begin{align}\boldsymbol{p}^{N{H}_{n}^{{\rm def}}}=\boldsymbol{p}^{NB}+\boldsymbol{p}^{B{R_{n}}}+\boldsymbol{p}^{{R_{n}}{H_{n}}}+{\boldsymbol{d}}_{n}^{{\rm rod}} \end{align}

where p ^{BR _n} is the position vector of the flexible rod’s mounting end relative to the platform’s center of mass, representing a constant vector denoted as p ^{BR _n} = $A_{n}^{\textrm{rod}}$ ${\boldsymbol{b}_{1}}$ + $B_{n}^{\textrm{rod}}$ ${\boldsymbol{b}_{2}}$ + $C_{n}^{\textrm{rod}}$ ${\boldsymbol{b}_{3}}$ (n = 1,⋯, N); p ^{R _n H _n} is the position vector of the mass element d $m_{n}^{\textrm{rod}}$ of the flexible rod in its undeformed state within the rod’s local coordinate frame, represented as p ^{R _n H _n} = l _n r _n1. Combining with Eq. (15), the coordinate matrix obtained in body-fixed coordinate frame B: { ${\boldsymbol{b}_{1}}$ , ${\boldsymbol{b}_{2}}$ , ${\boldsymbol{b}_{3}}$ } is:

(18) \begin{align}{}^{B}{\boldsymbol{p}}{_{}^{{R_{n}}{H_{n}}}}=\left({}^{B}{\boldsymbol{C}}{_{}^{R_n}}\right)^{\!\!\!{\rm T}}\left[l_{n},0,0\right]^{{\rm T}} \end{align}

As shown in Figure 6, $\boldsymbol{d}_{n}^{\textrm{rod}}$ represents the displacement of a mass element on a flexible rod. During displacement, this element undergoes lateral bending deformation, causing a projected contraction along the r _n1 direction. We define:

(19) \begin{align}{\boldsymbol{d}}_{n}^{{\rm rod}}=-{d}_{n1}^{{\rm rod}}\boldsymbol{r}_{n1}+{d}_{n2}^{{\rm rod}}\boldsymbol{r}_{n2}+{d}_{n3}^{{\rm rod}}\boldsymbol{r}_{n3} \end{align}

where ${d}_{n1}^{\textrm{rod}}$ represents the contraction displacement of the mass element along the r _n1 direction, which is a non-negative value; while ${d}_{n2}^{\textrm{rod}}$ and ${d}_{n3}^{\textrm{rod}}$ denote the bending deformation deflection of the mass element in the r _n2 and r _n3 directions, respectively. Using the assumed modal discretization method to extract the first three bending modes from yields the following expressions:

(20) \begin{align} \begin{cases} {d}_{n2}^{{\rm rod}}\left(l_{n},t\right)=\sum\limits_{j=1}^{3}\rho _{n2j}\left(t\right)\Psi _{j}\left(l_{n}\right)\\[5pt] {d}_{n3}^{{\rm rod}}\left(l_{n},t\right)=\sum\limits_{j=1}^{3}\rho _{n3j}\left(t\right)\Psi _{j}\left(l_{n}\right) \end{cases} \end{align}

where Ψ _j is the j _th mode shape function for a uniform cross-section cantilever beam; $\rho$ _n2j (t) and $\rho$ _n3j (t) are time-dependent modal coordinates. For convenience, denote:

\begin{align*}\rho _{n21}=\rho _{n1},\rho _{n22}=\rho _{n2},\rho _{n23}=\rho _{n3},\rho _{n31}=\rho _{n4},\rho _{n32}=\rho _{n5},\rho _{n33}=\rho _{n6} \end{align*}

The generalized coordinate vector for the flexible rod is defined as $\rho$ _n , and $\rho$ _n is expressed as:

(21) \begin{align}\boldsymbol{\rho }_{n}=\left[\rho _{n1},\rho _{n2},\rho _{n3},\rho _{n4},\rho _{n5},\rho _{n6}\right]^{{\rm T}}(n=1,\cdots , N) \end{align}

In the local coordinate system R _n : { r _n1, r _n2, r _n3}of the appendage rod, axial extensional deformation is neglected, and only the axial contraction caused by transverse bending deformation of the flexible rod is considered. Thus, for any mass element d $m_{n}^{\textrm{rod}}$ of the n _th flexible rod appendage (n = 1,⋯, N) in its undeformed state, its coordinates are [l _n ,0,0]^T. Based on the invariance of the arc length before and after deformation, the following equality relation is obtained:

(22) \begin{align}l_{n}=\int_{0}^{l_{n}-{d}_{n1}^{{\rm rod}}}\left[1+\left({\partial {d}_{n2}^{{\rm rod}}\left(\tau _{n},t\right) \over \partial \tau _{n}}\right)^{2}+\left({\partial {d}_{n3}^{{\rm rod}}\left(\tau _{n},t\right) \over \partial \tau _{n}}\right)^{2}\right]^{{1 \over 2}}{\rm d}\tau _{n},l_{n}\in \left[0,{L}_{n}^{{\rm rod}}\right] \end{align}

Under small transverse bending deformation, the following approximation can be adopted:

(23) \begin{align}\left[1+\left({\partial {d}_{n2}^{{\rm rod}}\left(\tau _{n},t\right) \over \partial \tau _{n}}\right)^{2}+\left({\partial {d}_{n3}^{{\rm rod}}\left(\tau _{n},t\right) \over \partial \tau _{n}}\right)^{2}\right]^{{1 \over 2}}\approx 1+{1 \over 2}\left[\left({\partial {d}_{n2}^{{\rm rod}}\left(\tau _{n},t\right) \over \partial \tau _{n}}\right)^{2}+\left({\partial {d}_{n3}^{{\rm rod}}\left(\tau _{n},t\right) \over \partial \tau _{n}}\right)^{2}\right] \end{align}

Substituting τ _n into Eq. (20) yields:

(24) \begin{align} \begin{cases} {d}_{n2}^{{\rm rod}}\left(\tau _{n},t\right)=\sum\limits_{j=1}^{3}\rho _{n2j}\left(t\right)\Psi _{j}\left(\tau _{n}\right)\\[8pt] {d}_{n3}^{{\rm rod}}\left(\tau _{n},t\right)=\sum\limits_{j=1}^{3}\rho _{n3j}\left(t\right)\Psi _{j}\left(\tau _{n}\right) \end{cases} \end{align}

Upon substituting Eqs. (23) and (24) into Eq. (22), one obtains:

(25) \begin{align}\begin{array}{c} {d}_{n1}^{{\rm rod}}=\dfrac{1 }{ 2}\int_{0}^{l_{n}-{d}_{n1}^{{\rm rod}}} \left[ \left( \dfrac{\partial \sum\limits_{j=1}^{3}\rho _{n2j}\left(t\right)\Psi _{j}\left(\tau _{n}\right) }{ \partial \tau _{n}} \right)^{2}+ \left( {\dfrac{\partial \sum\limits_{j=1}^{3}\rho _{n3j}\left(t\right)\Psi _{j}\left(\tau _{n}\right) }{ \partial \tau _{n}}}\right)^{2}\right]{\rm d}\tau _{n}\\[30pt] =\dfrac{1 }{ 2}\left({\boldsymbol{\rho }_{n2}}^{{\rm T}}\boldsymbol{G}_{2}\boldsymbol{\rho }_{n2}+{\boldsymbol{\rho }_{n3}}^{{\rm T}}\boldsymbol{G}_{3}\boldsymbol{\rho }_{n3}\right) \quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\quad\end{array} \end{align}

Upon assembling the matrices, the projected contraction of the geometrically nonlinear flexible cantilever beam admits the following compact representation:

(26) \begin{align}{d}_{n1}^{{\rm rod}}\left(l_{n},t\right)={1 \over 2}{\boldsymbol{\rho }_{n}}^{{\rm T}}\boldsymbol{G}_{n}\boldsymbol{\rho }_{n} \end{align}

where G _n is a 6 × 6 positive definite coefficient matrix, representing the geometric nonlinear projection contraction effect of the flexible rod.

Applying coordinate transformation via Eq. (15), the coordinate matrix of $\boldsymbol{d}_{n}^{\textrm{rod}}$ in the spacecraft body-fixed coordinate frame yields:

(27) \begin{align}{}^{B}{\boldsymbol{d}}{_{n}^{rod}}=\left({}^{B}{\boldsymbol{C}}{_{}^{Rn}}\right)^{{\rm T}}\left[-u_{n1},u_{n2},u_{n3}\right]^{{\rm T}} \end{align}

The generalized speed of the n _th(n = 1,⋯, N) flexible rod is now defined as:

(28) \begin{align}\boldsymbol{u}_{3k+6n+i}=\dot{\rho }_{ni}\left(i=1\sim 6,n=1,\cdots ,N\right) \end{align}

Finally, based on the relevant definitions, taking the derivative of $\boldsymbol{p}^{{NH}_{n}^{\textrm{def}}}$ with respect to time in the inertial reference frame N: { n ₁, n ₂, n ₃} yields the absolute velocity of any mass element d $m_{n}^{\textrm{rod}}$ on the flexible rod relative to the inertial reference frame:

(29) \begin{align}{}^{N}\boldsymbol{v}^{H_{n}^{\text{def}}} &= \frac{{}^{N}\text{d}}{\text{d}t} \boldsymbol{p}^{NH_{n}^{\text{def}}} = \frac{{}^{N}\text{d}}{\text{d}t} \left( \boldsymbol{p}^{NB} + \boldsymbol{p}^{BR_{n}} + \boldsymbol{p}^{R_{n}H_{n}} + \boldsymbol{d}_{n}^{\text{rod}} \right) \nonumber \\[3pt]&= {}^{N}\boldsymbol{v}^{B^{*}} + {}^{N}\boldsymbol{\omega}^{B} \times \boldsymbol{p}^{BR_{n}} + {}^{N}\boldsymbol{\omega}^{R_{n}} \times \left( \boldsymbol{p}^{R_{n}H_{n}} + \boldsymbol{d}_{n}^{\text{rod}} \right) + \frac{{}^{R_{n}}\text{d}}{\text{d}t} \boldsymbol{d}_{n}^{\text{rod}} \end{align}

2.3.4 Partial velocities and partial angular velocities

The liquid-filled spacecraft with flexible rod studied herein possesses 18 degrees of freedom. In the Kane method, the r _th partial velocity or partial angular velocity is defined as the coefficient of the r _th generalized speed in the corresponding velocity or angular velocity expression.

Spacecraft Rigid hub Partial Velocity and Partial Angular Velocities:

${}^{N}{\boldsymbol{v}_{1}^{\,B^{*}}}$ = ${\boldsymbol{b}_{1}}$ , ${}^{N}{\boldsymbol{v}}_{2}^{\,B^{*}}$ = ${\boldsymbol{b}_{2}}$ , ${}^{N}{\boldsymbol{v}}_{3}^{\,B^{*}}$ = ${\boldsymbol{b}_{3}}$

^N ${\boldsymbol{\omega}}_{1}^{\,B}$ = ${\boldsymbol{b}_{1}}$ , ^N ${\boldsymbol{\omega}}_{2}^{\,B}$ = ${\boldsymbol{b}_{2}}$ , ^N ${\boldsymbol{\omega}}_{3}^{\,B}$ = ${\boldsymbol{b}_{3}}$

Partial Velocity in the Liquid Sloshing Equivalent Center of Mass Constraint Model:

${}^{N}\boldsymbol{v}_{1}^{\,P_{1}}$ = ${\boldsymbol{b}_{1}}$ , ${}^{N}\boldsymbol{v}_{2}^{\,P_{1}}$ = ${\boldsymbol{b}_{2}}$ , ${}^{N}\boldsymbol{v}_{3}^{\,P_{1}}$ = ${\boldsymbol{b}_{3}}$

\begin{align*}{}^{N}{\boldsymbol{v}}{_{4}^{\,P_{1}}}=\boldsymbol{b}_{3}(Y_{1}+y_{1})-\boldsymbol{b}_{2}(Z_{1}+z_{1}) \\[-23pt] \end{align*}

\begin{align*}{}^{N}{\boldsymbol{v}}{_{5}^{\,P_{1}}}=\boldsymbol{b}_{1}(Z_{1}+z_{1})-\boldsymbol{b}_{3}(X_{1}+x_{1}) \\[-23pt] \end{align*}

\begin{align*}{}^{N}{\boldsymbol{v}}{_{6}^{\,P_{1}}}=\boldsymbol{b}_{2}(X_{1}+x_{1})-\boldsymbol{b}_{1}(Y_{1}+y_{1}) \end{align*}

${}^{N}\boldsymbol{v}_{7}^{\,P_{1}}$ = ${\boldsymbol{b}_{1}}$ , ${}^{N}\boldsymbol{v}_{8}^{\,P_{1}}$ = ${\boldsymbol{b}_{2}}$ , ${}^{N}\boldsymbol{v}_{9}^{\,P_{1}}$ = ${\boldsymbol{b}_{3}}$

${}^{N}\boldsymbol{v}_{1}^{\,P_{2}}$ = ${\boldsymbol{b}_{1}}$ , ${}^{N}\boldsymbol{v}_{2}^{\,P_{2}}$ = ${\boldsymbol{b}_{2}}$ , ${}^{N}\boldsymbol{v}_{3}^{\,P_{2}}$ = ${\boldsymbol{b}_{3}}$

\begin{align*}{}^{N}{\boldsymbol{v}}{_{4}^{\,P_{2}}}=\boldsymbol{b}_{3}(Y_{2}+y_{2})-\boldsymbol{b}_{2}(Z_{2}+z_{2}) \\[-23pt] \end{align*}

\begin{align*}{}^{N}{\boldsymbol{v}}{_{5}^{\,P_{2}}}=\boldsymbol{b}_{1}(Z_{2}+z_{2})-\boldsymbol{b}_{3}(X_{2}+x_{2}) \\[-23pt] \end{align*}

\begin{align*}{}^{N}{\boldsymbol{v}}{_{6}^{\,P_{2}}}=\boldsymbol{b}_{2}(X_{2}+x_{2})-\boldsymbol{b}_{1}(Y_{2}+y_{2}) \end{align*}

${}^{N}\boldsymbol{v}_{10}^{\,P_{1}}$ = ${\boldsymbol{b}_{1}}$ , ${}^{N}\boldsymbol{v}_{11}^{\,P_{1}}$ = ${\boldsymbol{b}_{2}}$ , ${}^{N}\boldsymbol{v}_{12}^{\,P_{1}}$ = ${\boldsymbol{b}_{3}}$

Flexible Rod Displacement Partial Velocity and Partial Angular Velocities:

Rigid hub Rotation of flexible rod:

^N ${\boldsymbol{\omega}}_{4}^{\,R_{1}}$ = ${\boldsymbol{b}_{1}}$ , ^N ${\boldsymbol{\omega}}_{5}^{\,R_{1}}$ = ${\boldsymbol{b}_{2}}$ , ^N ${\boldsymbol{\omega}}_{6}^{\,R_{1}}$ = ${\boldsymbol{b}_{3}}$ .

Flexible rod translation:

${}^{N}{\boldsymbol{v}}_{1}^{\,H^{\textrm{def}}_{1}}$ = ${\boldsymbol{b}_{1}}$ , ${}^{N}{\boldsymbol{v}}_{2}^{\,H^{\textrm{def}}_{1}}$ = ${\boldsymbol{b}_{2}}$ , ${}^{N}{\boldsymbol{v}}_{3}^{\,H^{\textrm{def}}_{1}}$ = ${\boldsymbol{b}_{3}}$ .

Vibration of flexible rod:

\begin{align*}\begin{array}{l} {}^{N}{\boldsymbol{v}}{_{4}^{H_{1}^{{\rm def}}}}=\left[\left(l_{1}-\dfrac{1}{2}\sum\limits_{i=1}^{6}\sum\limits_{j=1}^{6}G_{1ij}\rho _{1i}\rho _{1j}\right)\sin \alpha _{1} \cos \beta _{1}-\left(\sum\limits_{i=1}^{3}\rho _{1i}\Psi _{i}\right)\sin \alpha _{1}\sin \beta _{1}-\left(\sum\limits_{i=1}^{3}{{\rho}_{1(3+i)}}\Psi _{i}\right)\cos \alpha _{1}\right]\boldsymbol{b}_{2}\\[15pt] \quad \quad \quad \ +\left[{B}_{1}^{{\rm rod}}+\left(l_{1}-\dfrac{1}{2}\sum\limits_{i=1}^{6}\sum\limits_{j=1}^{6}G_{1ij}\rho _{1i}\rho _{1j}\right)\sin \beta _{1}+\left(\sum\limits_{i=1}^{3}\rho _{1i}\Psi _{i}\right)\cos \beta _{1}\right]\boldsymbol{b}_{3} \end{array} \end{align*}

\begin{align*}\begin{array}{l} {}^{N}{\boldsymbol{v}}{_{5}^{H_{1}^{{\rm def}}}}=-\left[\left(l_{1}-\dfrac{1}{2}\sum\limits_{i=1}^{6}\sum\limits_{j=1}^{6}G_{1ij}\rho _{1i}\rho _{1j}\right)\sin \alpha _{1}\cos \beta _{1}-\left(\sum\limits_{i=1}^{3}\rho _{1i}\Psi _{i}\right)\sin \alpha _{1}\sin \beta _{1}-\left(\sum\limits_{i=1}^{3}\rho _{1(3+i)}\Psi _{i}\right)\cos \alpha _{1}\right]\boldsymbol{b}_{1}\\[15pt] \quad \quad \quad \ -\left[{A}_{1}^{{\rm rod}}+\left(l_{1}-\dfrac{1}{2}\sum\limits_{i=1}^{6}\sum\limits_{j=1}^{6}G_{1ij}\rho _{1i}\rho _{1j}\right)\cos \alpha _{1}\cos \beta _{1}-\left(\sum\limits_{i=1}^{3}\rho _{1i}\Psi _{i}\right)\cos \alpha _{1}\sin \beta _{1}+\left(\sum\limits_{i=1}^{3}\rho _{1(3+i)}\Psi _{i}\right)\sin \alpha _{1}\right]\boldsymbol{b}_{3} \end{array} \end{align*}

\begin{align*}\begin{array}{l} {}^{N}{\boldsymbol{v}}{_{6}^{H_{1}^{{\rm def}}}}=-\left[{B}_{1}^{{\rm rod}}+\left(l_{1}-\dfrac{1}{2}\sum\limits_{i=1}^{6}\sum\limits_{j=1}^{6}G_{1ij}\rho _{1i}\rho _{1j}\right)\sin \beta _{1}+\left(\sum\limits_{i=1}^{3}\rho _{1i}\Psi _{i}\right)\cos \beta _{1}\right]\boldsymbol{b}_{1}\\[15pt]\quad \quad\quad \ +\left[{A}_{1}^{{\rm rod}}+\left(l_{1}-\dfrac{1}{2}\sum\limits_{i=1}^{6}\sum\limits_{j=1}^{6}G_{1ij}\rho _{1i}\rho _{1j}\right)\cos \alpha _{1}\cos \beta _{1}-\left(\sum\limits_{i=1}^{3}\rho _{1i}\Psi _{i}\right)\cos \alpha _{1}\sin \beta _{1}+\left(\sum\limits_{i=1}^{3}\rho _{1(3+i)}\Psi _{i}\right)\sin \alpha _{1}\right]\boldsymbol{b}_{2} \end{array} \\[-23pt] \end{align*}

\begin{align*}\begin{array}{l} {}^{N}{\boldsymbol{v}}{_{13}^{H_{1}^{{\rm def}}}}=\left[-\left(\sum\limits_{j=1}^{6}G_{11j}\rho _{1j}\right)\cos \alpha _{1}\cos \beta _{1}-\Psi _{1}\cos \alpha _{1}\sin \beta _{1}\right]\boldsymbol{b}_{1}+\left[-\left(\sum\limits_{j=1}^{6}G_{11j}\rho _{1j}\right)\sin \beta _{1}+\Psi _{1}\cos \beta _{1}\right]\boldsymbol{b}_{2}\\[10pt] \quad \quad \quad \ +\left[\left(\sum\limits_{j=1}^{6}G_{11j}\rho _{1j}\right)\sin \alpha _{1}\cos \beta _{1}+\Psi _{1}\sin \alpha _{1}\sin \beta _{1}\right]\boldsymbol{b}_{3} \end{array} \\[-23pt] \end{align*}

\begin{align*}\begin{array}{l} {}^{N}{\boldsymbol{v}}{_{14}^{H_{1}^{{\rm def}}}}=\left[-\left(\sum\limits_{j=1}^{6}G_{12j}\rho _{1j}\right)\cos \alpha _{1}\cos \beta _{1}-\Psi _{2}\cos \alpha _{1}\sin \beta _{1}\right]\boldsymbol{b}_{1}+\left[-\left(\sum\limits_{j=1}^{6}G_{12j}\rho _{1j}\right)\sin \beta _{1}+\Psi _{2}\cos \beta _{1}\right]\boldsymbol{b}_{2}\\[15pt] \quad \quad \quad \ +\left[\left(\sum\limits_{j=1}^{6}G_{12j}\rho _{1j}\right)\sin \alpha _{1}\cos \beta _{1}+\Psi _{2}\sin \alpha _{1}\sin \beta _{1}\right]\boldsymbol{b}_{3} \end{array} \\[-23pt] \end{align*}

\begin{align*}\begin{array}{l} {}^{N}{\boldsymbol{v}}{_{15}^{H_{1}^{{\rm def}}}}=\left[-\left(\sum\limits_{j=1}^{6}G_{13j}\rho _{1j}\right)\cos \alpha _{1}\cos \beta _{1}-\Psi _{3}\cos \alpha _{1}\sin \beta _{1}\right]\boldsymbol{b}_{1}+\left[-\left(\sum\limits_{j=1}^{6}G_{13j}\rho _{1j}\right)\sin \beta _{1}+\Psi _{3}\cos \beta _{1}\right]\boldsymbol{b}_{2}\\[15pt] \quad \quad \quad \ +\left[\left(\sum\limits_{j=1}^{6}G_{13j}\rho _{1j}\right)\sin \alpha _{1}\cos \beta _{1}+\Psi _{3}\sin \alpha _{1}\sin \beta _{1}\right]\boldsymbol{b}_{3} \end{array} \\[-23pt] \end{align*}

\begin{align*}\begin{array}{l} {}^{N}{\boldsymbol{v}}{_{16}^{H_{1}^{{\rm def}}}}=\left[-\bigg(\sum\limits_{j=1}^{6}G_{14j}\rho _{1j}\bigg)\cos \alpha _{1}\cos \beta _{1}+\Psi _{1}\sin \alpha _{1}\right]\boldsymbol{b}_{1}+\left[-\bigg(\sum\limits_{j=1}^{6}G_{14j}\rho _{1j}\bigg)\sin \beta _{1}\right]\boldsymbol{b}_{2}\\[15pt] \quad \quad \quad \ +\left[\bigg(\sum\limits_{j=1}^{6}G_{14j}\rho _{1j}\bigg)\sin \alpha _{1}\cos \beta _{1}+\Psi _{1}\cos \alpha _{1}\right]\boldsymbol{b}_{3} \end{array}\\[-23pt] \end{align*}

\begin{align*}\begin{array}{l} {}^{N}{\boldsymbol{v}}{_{17}^{H_{1}^{{\rm def}}}}=\left[-\bigg(\sum\limits_{j=1}^{6}G_{15j}\rho _{1j}\bigg)\cos \alpha _{1}\cos \beta _{1}+\Psi _{2}\sin \alpha _{1}\right]\boldsymbol{b}_{1}+\left[-\bigg(\sum\limits_{j=1}^{6}G_{15j}\rho _{1j}\bigg)\sin \beta _{1}\right]\boldsymbol{b}_{2}\\[15pt] \quad \quad \quad \ +\left[\bigg(\sum\limits_{j=1}^{6}G_{15j}\rho _{1j}\bigg)\sin \alpha _{1}\cos \beta _{1}+\Psi _{2}\cos \alpha _{1}\right]\boldsymbol{b}_{3} \end{array}\\[-23pt] \end{align*}

\begin{align*}\begin{array}{l} {}^{N}{\boldsymbol{v}}{_{18}^{H_{1}^{{\rm def}}}}=\left[-\bigg(\sum\limits_{j=1}^{6}G_{16j}\rho _{1j}\bigg)\cos \alpha _{1}\cos \beta _{1}+\Psi _{3}\sin \alpha _{1}\right]\boldsymbol{b}_{1}+\left[-\bigg(\sum\limits_{j=1}^{6}G_{16j}\rho _{1j}\bigg)\sin \beta _{1}\right]\boldsymbol{b}_{2}\\[15pt] \quad \quad \quad \ +\left[\bigg(\sum\limits_{j=1}^{6}G_{16j}\rho _{1j}\bigg)\sin \alpha _{1}\cos \beta _{1}+\Psi _{3}\cos \alpha _{1}\right]\boldsymbol{b}_{3} \end{array} \end{align*}

2.3.5 Generalized inertial forces

By differentiating the corresponding velocities and angular velocities of each part of the spacecraft with respect to time, the corresponding accelerations and angular accelerations can be obtained. The absolute acceleration of the main rigid hub’s center of mass B ^*during translational motion and the absolute angular acceleration of the rigid hub during rotation, both relative to the inertial frame, are respectively given by Eqs. (30) and (31).

(30) \begin{align} \begin{array}{l} {}^{N}{\boldsymbol{a}}{_{}^{B^{*}}} = \dfrac{{}^{N}{{\rm d}}{}}{{\rm d}t}{}^{N}{\boldsymbol{v}}{_{}^{B^{*}}}={}^{N}{\boldsymbol{\omega }}{_{}^{B}}\times {}^{N}{\boldsymbol{v}}{_{}^{B^{*}}}+\dfrac{{}^{B}{{\rm d}}{} }{ {\rm d}t}{}^{N}{\boldsymbol{v}}{_{}^{B^{*}}}\\[8pt] \quad \quad =\boldsymbol{b}_{1}\left(\, \dot{\kern-1pt u}_{1}+u_{3}u_{5}-u_{2}u_{6}\right)\\[3pt] \quad \quad\quad +\boldsymbol{b}_{2}\left(\, \dot{\kern-1pt u}_{2}+u_{1}u_{6}-u_{3}u_{4}\right)\\[3pt] \quad \quad \quad +\boldsymbol{b}_{3}\left(\, \dot{\kern-1pt u}_{3}+u_{2}u_{4}-u_{1}u_{5}\right) \end{array} \\[-6pt] \nonumber \end{align}

(31) \begin{align} {}^{N}{\boldsymbol{\alpha }}{_{}^{B}}={{}^{N}{{\rm d}}{} \over {\rm d}t}{}^{N}{\boldsymbol{\omega }}{_{}^{B}}=\, \dot{\kern-1pt u}_{3+i}\cdot \boldsymbol{b}_{i}=\, \dot{\kern-1pt u}_{4}\boldsymbol{b}_{1}+\, \dot{\kern-1pt u}_{5}\boldsymbol{b}_{2}+\, \dot{\kern-1pt u}_{6}\boldsymbol{b}_{3} \\[12pt] \nonumber \end{align}

Equation (32) represents the absolute velocity of the liquid-filled tank, while the absolute acceleration and absolute angular velocity of the flexible rod are given by Eqs. (33) and (34), respectively.

(32) \begin{align} \begin{array}{l} {}^{N}{\boldsymbol{a}}{_{}^{P_{k}}}= \dfrac{{}^{N}{{\rm d}}{} }{ {\rm d}t}{}^{N}{\boldsymbol{v}}{_{}^{P_{k}}}={}^{N}{\boldsymbol{a}}{_{}^{B^{*}}}+{}^{N}{\boldsymbol{\omega }}{_{}^{B}}\times \left[{}^{N}{\boldsymbol{\omega }}{_{}^{B}}\times \left(\boldsymbol{p}^{{B^{*}}{T_{k}}}+\boldsymbol{p}^{{T_{k}}{P_{k}}}\right)+2 \dfrac{{}^{B}{{\rm d}}{} }{ {\rm d}t}\boldsymbol{p}^{{T_{k}}{P_{k}}}\right]+{}^{N}{\boldsymbol{\alpha }}{_{}^{B}}\times \left(\boldsymbol{p}^{{B^{*}}{T_{k}}}+\boldsymbol{p}^{{T_{k}}{P_{k}}}\right)+ \dfrac{{}^{B}{{\rm d}}{_{}^{2}} }{ {\rm d}t^{2}}\boldsymbol{p}^{{T_{k}}{P_{k}}}\\[8pt]\quad \quad =\boldsymbol{b}_{1}\left\{\, \dot{\kern-1pt u}_{1}+\left(Z_{k}+z_{k}\right)\, \dot{\kern-1pt u}_{5}-\left(Y_{k}+y_{k}\right)\, \dot{\kern-1pt u}_{6}+\, \dot{\kern-1pt u}_{7+3k}+u_{3}u_{5}-u_{2}u_{6}-\left(X_{k}+x_{k}\right){u_{5}}^{2}\right.\\[5pt]\quad \quad \quad \left.-\left(X_{k}+x_{k}\right){u_{6}}^{2}+\left(Y_{k}+y_{k}\right)u_{4}u_{5}+\left(Z_{k}+z_{k}\right)u_{4}u_{6}+2u_{5}u_{9+3k}-2u_{6}u_{8+3k}\right\}\\[5pt]\quad \quad \quad +\boldsymbol{b}_{2}\left\{\, \dot{\kern-1pt u}_{2}-\left(Z_{k}+z_{k}\right)\, \dot{\kern-1pt u}_{4}+\left(X_{k}+x_{k}\right)\, \dot{\kern-1pt u}_{6}+\, \dot{\kern-1pt u}_{8+3k}-u_{3}u_{4}+u_{1}u_{6}-\left(Y_{k}+y_{k}\right){u_{4}}^{2}\right.\\[5pt]\quad \quad \quad \left.-\left(Y_{k}+y_{k}\right){u_{6}}^{2}+\left(X_{k}+x_{k}\right)u_{4}u_{5}+\left(Z_{k}+z_{k}\right)u_{5}u_{6}-2u_{4}u_{9+3k}+2u_{6}u_{7+3k}\right\}\\[5pt]\quad \quad \quad +\boldsymbol{b}_{3}\left\{\, \dot{\kern-1pt u}_{3}+\left(Y_{k}+y_{k}\right)\, \dot{\kern-1pt u}_{4}-\left(X_{k}+x_{k}\right)\, \dot{\kern-1pt u}_{5}+\, \dot{\kern-1pt u}_{9+3k}+u_{2}u_{4}-u_{1}u_{5}-\left(Z_{k}+z_{k}\right){u_{4}}^{2}\right.\\[5pt]\quad \quad \quad \left.-\left(Z_{k}+z_{k}\right){u_{5}}^{2}+\left(X_{k}+x_{k}\right)u_{4}u_{6}+\left(Y_{k}+y_{k}\right)u_{5}u_{6}+2u_{4}u_{8+3k}-2u_{5}u_{7+3k}\right\}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace}{\enspace} \end{array} \nonumber \\[-20pt] \end{align}

(33) \begin{align} \begin{array}{l} {}^{N}{\boldsymbol{a}}{_{}^{H_{1}^{{\rm def}}}}= \dfrac{{}^{N}{{\rm d}}{} }{{\rm d}t}{}^{N}{\boldsymbol{v}}{_{}^{H_{1}^{{\rm def}}}}\\[10pt] \quad \quad \ = \dfrac{{}^{N}{{\rm d}}{} }{ {\rm d}t}\left[{}^{N}{\boldsymbol{v}}{_{}^{B^{*}}}+{}^{N}{\boldsymbol{\omega }}{_{}^{B}}\times \boldsymbol{p}^{B{R_{1}}}+{}^{N}{\boldsymbol{\omega }}{^{R_{1}}}\times \left(\boldsymbol{p}^{{R_{1}}{H_{1}}}+{\boldsymbol{d}}_{1}^{{\rm rod}}\right)+ \dfrac{{}^{R_{1}}{{\rm d}}{} }{{\rm d}t}{\boldsymbol{d}}_{1}^{{\rm rod}}\right]\\[12pt] \quad \quad \ ={}^{N}{\boldsymbol{a}}{^{B^{*}}}+{}^{N}{\boldsymbol{\alpha }}{_{}^{B}}\times \boldsymbol{p}^{B{R_{1}}}+{}^{N}{\boldsymbol{\omega }}{_{}^{B}}\times \left({}^{N}{\boldsymbol{\omega }}{_{}^{B}}\times \boldsymbol{p}^{B{R_{1}}}\right)+{}^{N}{\boldsymbol{\alpha }}{_{}^{R_{1}}}\times \left(\boldsymbol{p}^{{R_{1}}{H_{1}}}+{\boldsymbol{d}}_{1}^{{\rm rod}}\right)\\[8pt] \quad \quad \quad +{}^{N}{\boldsymbol{\omega }}{^{R_{1}}}\times \left[{}^{N}{\boldsymbol{\omega }}{_{}^{R_{1}}}\times \left(\boldsymbol{p}^{{R_{1}}{H_{1}}}+{\boldsymbol{d}}_{1}^{{\rm rod}}\right)+2\dfrac{{}^{R_{1}}{{\rm d}}{}}{{\rm d}t}{\boldsymbol{d}}_{1}^{{\rm rod}}\right]+ \dfrac{{}^{R_{1}}{{\rm d}}{_{}^{2}}}{{\rm d}t^{2}}{\boldsymbol{d}}_{1}^{{\rm rod}} \end{array} \\[-20pt] \nonumber \end{align}

(34) \begin{align} \begin{array}{ll} {}^{N}{\boldsymbol{\alpha }}{_{}^{R_{1}}} & =\dfrac{{}^{N}{{\rm d}}{}}{ {\rm d}t}\left({}^{N}{\boldsymbol{\omega }}{_{}^{B}}+\dot{\alpha }_{1}{\boldsymbol{r}}_{12}^{0}+\dot{\beta }_{1}\boldsymbol{r}_{13}\right)\\[3pt] & ={}^{N}{\boldsymbol{\alpha }}{_{}^{B}}+{}^{N}{\boldsymbol{\omega }}{_{}^{R_{1}^{0}}}\times \left(\dot{\alpha }_{1}{\boldsymbol{r}}_{12}^{0}\right)+ \dfrac{{}^{R_{1}^0}{{\rm d}}{}}{{\rm d}t}\left(\dot{\alpha }_{1}{\boldsymbol{r}}_{12}^{0}\right)+{}^{N}{\boldsymbol{\omega }}{^{R_{1}}}\times \left(\dot{\beta }_{1}\boldsymbol{r}_{13}\right)+ \dfrac{{}^{R_{1}}{{\rm d}}{}}{{\rm d}t}\left(\dot{\beta }_{1}\boldsymbol{r}_{13}\right) \end{array} \\[3pt] \nonumber \end{align}

The inertial forces acting on the spacecraft discussed herein are:

(35) \begin{align}\begin{array}{l} {F}_{r}^{*}=-m_{B}{}^{N}{\boldsymbol{a}}{^{B^{*}}}\cdot {}^{N}{\boldsymbol{v}}{_{r}^{B^{*}}}-\sum\limits_{k=1}^{2}{m}_{k}^{{\rm liq}} \, {}^{N}{\boldsymbol{a}}{^{P_{k}}}\cdot {}^{N}{\boldsymbol{v}}{_{r}^{P_{k}}}-\left(\boldsymbol{I}^{B}\cdot {}^{N}{\boldsymbol{\alpha }}{_{}^{B}}+{}^{N}{\boldsymbol{\omega }}{_{}^{B}}\times \boldsymbol{I}^{B}\cdot {}^{N}{\boldsymbol{\omega }}{_{}^{B}}\right)\cdot {}^{N}{\boldsymbol{\omega }}{_{r}^{B}}\\[9pt] \quad \quad \, - \left(\boldsymbol{I}^{{\rm rod}}\cdot {}^{N}{\boldsymbol{\alpha }}{_{}^{R_{1}}}+{}^{N}{\boldsymbol{\omega }}{_{}^{R_{1}}}\times \boldsymbol{I}^{{\rm rod}}\cdot {}^{N}{\boldsymbol{\omega }}{^{R_{1}}}\right)\cdot {}^{N}{\boldsymbol{\omega }}{_{r}^{R_{1}}}-\int _{{m}_{1}^{{\rm rod}}}{}^{N}{\boldsymbol{a}}{_{}^{{H_{1}}^{{\rm def}}}}\cdot {}^{N}{\boldsymbol{v}_{r}^{{H}_{1}^{{\rm def}}}}{\rm d}{m}_{1}^{{\rm rod}} \end{array} \end{align}

where m _B is the mass of the spacecraft’s main rigid hub; ${m}_{k}^{\textrm{liq}}$ is the mass of liquid propellant in the k _th tank; I ^B is the second-order inertia tensor of the spacecraft’s main rigid hub relative to the body-fixed coordinate frame. I ^rod is the second-order inertia tensor of the flexible rod relative to its local coordinate frame; for computational simplicity, its deformation effects are neglected in this paper; ${m}_{1}^{\textrm{rod}}$ is the mass of the rod.

2.3.6 Generalized force

Consider the generalized force ${F}_{r}$ exerted on the spacecraft in low-gravity or microgravity environments, arising from gravitational forces and the constraint force from the equivalent liquid center-of-mass constraint surface model within the tank. Its calculation formula is:

(36) \begin{align}F_{r}={\boldsymbol{F}}_{g}^{B^{*}}\cdot {}^{N}{\boldsymbol{v}}{_{r}^{B^{*}}}+\sum\limits_{k=1}^{2}\left({\boldsymbol{F}}_{g}^{P_{k}}+{\boldsymbol{F}_{c}}^{{P_{{\rm k}}}}\right)\cdot {}^{N}{\boldsymbol{v}}{_{r}^{P_{k}}}+\boldsymbol{F}_{c}\cdot {}^{N}{\boldsymbol{v}}{_{r}^{B^{*}}}+\boldsymbol{T}_{c}\cdot {}^{N}{\boldsymbol{\omega }}{_{r}^{B}}+\int _{{m}_{1}^{{\rm rod}}}\left({\boldsymbol{F}}_{g}^{{H}_{1}^{{\rm def}}}\cdot {}^{N}{\boldsymbol{v}}{_{r}^{H_{1}^{{\rm def}}}}\right)+{F}_{12+i}^{{\rm rod}} \end{align}

where $\boldsymbol{F}_{g}^{B^{*}}$ represents the gravitational force acting on the spacecraft’s main rigid hub center of mass, expressed as Eq. (37); $\boldsymbol{F}_{c}{\enspace}$ the control thrust balancing the gravitational force acting on the spacecraft, expressed as Eq. (38); T _c represents its attitude control torque, expressed as Eq. (39); $\boldsymbol{F}_{g}^{P_{k}}$ the gravitational force acting on the equivalent center of mass of the liquid in the k _th tank is expressed as Eq. (40).

(37) \begin{align} {\boldsymbol{F}}_{g}^{B^{*}}=-m_{B}g\boldsymbol{n}_{1}=-m_{B}g\sum\limits_{j=1}^{3}{}^{N}{C}{_{j1}}{{}^{B}\boldsymbol{b}_{j}} \\[-10pt] \nonumber \end{align}

(38) \begin{align} \boldsymbol{F}_{c}=F_{{c_{1}}}\boldsymbol{b}_{1}+F_{{c_{2}}}\boldsymbol{b}_{2}+F_{{c_{3}}}\boldsymbol{b}_{3}\quad \quad \\[-5pt] \nonumber \end{align}

(39) \begin{align} \boldsymbol{T}_{c}=T_{{c_{1}}}\boldsymbol{b}_{1}+T_{{c_{2}}}\boldsymbol{b}_{2}+T_{{c_{3}}}\boldsymbol{b}_{3} \quad \quad \\[-9pt] \nonumber \end{align}

(40) \begin{align} {\boldsymbol{F}}_{g}^{P_{k}}=-{m}_{k}^{liq}g\boldsymbol{n}_{1}=-{m}_{k}^{liq}g\sum\limits_{j=1}^{3}{}^{N}{C}{_{j1}}{{}^{B}\boldsymbol{b}_{j}} \\[10pt] \nonumber \end{align}

where $\boldsymbol{F}_{c}^{\,P_{k}}$ represents the constraint force exerted on the equivalent center of mass within the tank by the corresponding center-of-mass constraint surface. When the equivalent center of mass moves within the center-of-mass constraint surface, the constraint equation ensuring its motion on the ellipsoid is:

(41) \begin{align}\left(\sum\limits_{i=1}^{3}u_{3+3k+i}\boldsymbol{b}_{i}\right)\cdot \hat{\boldsymbol{n}}_{k}=0 \end{align}

where $\hat{\boldsymbol{n}}_{k}(k=1,2)$ is the outward normal vector of the equivalent center of mass moving along the center-of-mass constraint surface inside the tank. Taking the total differential of Eq. (1) yields:

(42) \begin{align}\hat{\boldsymbol{n}}_{k}={1 \over 2}\left(\boldsymbol{b}_{1}\left.{\partial S_{k} \over \partial x}\right| _{x={x_{k}}}+\boldsymbol{b}_{2}\left.{\partial S_{k} \over \partial y}\right| _{y={y_{k}}}+\boldsymbol{b}_{3}\left.{\partial S_{k} \over \partial z}\right| _{z={z_{k}}}\right)={x_{k} \over {a_{k}}^{2}}\boldsymbol{b}_{1}+{y_{k} \over {b_{k}}^{2}}\boldsymbol{b}_{2}+{z_{k} \over {b_{k}}^{2}}\boldsymbol{b}_{3} \end{align}

Differentiating Eq. (41) with respect to time yields the constraint equation in acceleration form, which is presented as Eq. (43). It is worth clarifying that the above constraint is holonomic. This constraint is directly derived from Eq. (1), which defines the requirement that the mass center must move on the ellipsoidal surface. In other words, the acceleration-form constraint equation presented above is merely the second-order time differentiation of the original geometric constraint.

(43) \begin{align}\left(\sum\limits_{i=1}^{3}\, \dot{\kern-1pt u}_{3+3k+i}\boldsymbol{b}_{i}\right)\cdot \hat{\boldsymbol{n}}_{k}+\left(\sum\limits_{i=1}^{3}u_{3+3k+i}\boldsymbol{b}_{i}\right)\cdot {{}^{B}{{\rm d}}{} \over {\rm d}t}\hat{\boldsymbol{n}}_{k}=0 \end{align}

It can be seen that the above equation is a first-order linear expression in terms of the generalized speed, which can be represented in the following matrix form:

(44) \begin{align}\boldsymbol{A}\, \dot{\boldsymbol{\kern-1pt U}}+\boldsymbol{B}={\bf 0} \end{align}

where U = [u ₁,⋯,u ₁₈]^T is the system state variable matrix of dimensions 18 × 1, and simultaneously $\, \dot{\boldsymbol{\kern-1pt U}}=[\, \dot{\kern-1pt u}_{1},\cdots ,\, \dot{\kern-1pt u}_{18}]^{{\rm T}}$ .

The constraint force of the center constraint surface on the equivalent centroid point can be decomposed into normal support force and tangential friction force, namely:

(45) \begin{align}{\boldsymbol{F}}_{c}^{P_{k}}=-\mu _{k}\left(\sum\limits_{i=1}^{3}u_{3+3k+i}\boldsymbol{b}_{i}\right)-\lambda _{k}{\hat{\boldsymbol{n}}_{k} \over \left| \hat{\boldsymbol{n}}_{k}\right| }\left(k=1\sim 2\right) \end{align}

where $\mu _{k}{\enspace}(k=1\sim 2)$ is the equivalent friction characteristic parameter of the liquid, characterizing the dissipative properties of the equivalent center-of-mass point motion in the liquid sloshing equivalent center-of-mass constraint surface model. Its value is influenced by the coupling effects of external excitation (e.g., gravitational environment), tank structural parameters (tank dimensions), liquid fill ratio, and liquid physical properties (such as liquid density and viscosity). It can be quantified through liquid sloshing experiments or CFD numerical simulations. $\lambda _{k}{\enspace}(k=1\sim 2)$ is the Lagrangian undetermined multiplier, representing the scalar value of the normal support force acting on the equivalent center of mass point.

$\boldsymbol{F}_{g}^{H_{1}^{\textrm{def}}}$ The specific expression is:

(46) \begin{align}{\boldsymbol{F}}_{g}^{{H}_{1}^{{\rm def}}}=-\boldsymbol{n}_{1}g{\rm d}{m}_{1}^{{\rm rod}}=-\left(\sum\limits_{j=1}^{3}{}^{N}{C}{_{j1}^{B}}\boldsymbol{b}_{j}\right)g{\rm d}{m}_{1}^{{\rm rod}} \end{align}

${{F}_{12+i}^{\textrm{rod}}}$ are the GAFs due to the material elasticity of the flexible rod, which can be calculated by:

(47) \begin{align}{F}_{12+i}^{{\rm rod}}=-\sum\limits_{j=1}^{6}{K}_{ij}^{{\rm rod}}\boldsymbol{\rho }_{ij}\left(i=1\sim 6\right) \end{align}

where $K_{ij}^{\textrm{rod}}$ is the element of the modal stiffness matrix ${\boldsymbol{K}_{1}^{\textrm{rod}}}$ . The modal stiffness matrix ${\boldsymbol{K}_{1}^{\textrm{rod}}}$ is defined by calculating the strain energy stored during lateral vibration of the flexible rod, i.e.,

(48) \begin{align}\begin{array}{l} {V}_{1}^{{\rm rod}}=\dfrac{1 }{ 2}\displaystyle\int_{0}^{{L}_{1}^{{\rm rod}}}\left[{E}_{1}^{{\rm rod}}{I}_{13}^{{\rm rod}}\left( \dfrac{\partial ^{2}{d}_{12}^{{\rm rod}}(x,t) }{ \partial x^{2}}\right)^{2}+{E}_{1}^{{\rm rod}}{I}_{12}^{{\rm rod}}\left( \dfrac{\partial ^{2}{d}_{13}^{{\rm rod}}(x,t) }{\partial x^{2}}\right)^{2}\right]{\rm d}x\\[12pt] \quad \quad =\dfrac{1 }{ 2}{\boldsymbol{\rho }}_{1}^{{\rm T}}{\boldsymbol{K}}_{1}^{{\rm rod}}\boldsymbol{\rho }_{1} \end{array} \end{align}

where E ₁ ^rod is the Young’s modulus of the flexible rod, I ₁₂ ^rod and I ₁₃ ^rod are the product of the moments of inertia of the flexible rod along the r ₁₂ and r ₁₃ axes, respectively, and ${\boldsymbol{K}_{1}^{\textrm{rod}}}$ is the 6 × 6 modal stiffness matrix:

\begin{align*}{\boldsymbol{K}}_{1}^{{\rm rod}}=\left[\begin{array}{ll} {\boldsymbol{K}}_{12}^{{\rm rod}} &\quad \boldsymbol{0}_{3\times 3}\\[3pt] \boldsymbol{0}_{3\times 3} &\quad {\boldsymbol{K}}_{13}^{{\rm rod}} \end{array}\right] \\[-10pt] \end{align*}

where $\begin{cases} {\boldsymbol{K}}_{12}^{{\rm rod}}\left(i,j\right)={E}_{1}^{{\rm rod}}{I}_{13}^{{\rm rod}}\displaystyle\int_{0}^{{L}_{1}^{{\rm rod}}} \dfrac{\partial ^{2}\Psi _{i} }{ \partial x^{2}} \dfrac{\partial ^{2}\Psi _{j}}{\partial x^{2}}{\rm d}x\\[8pt] {\boldsymbol{K}}_{13}^{{\rm rod}}\left(i,j\right)={E}_{1}^{{\rm rod}}{I}_{12}^{{\rm rod}}\displaystyle\int_{0}^{{L}_{1}^{{\rm rod}}} \dfrac{\partial ^{2}\Psi _{i} }{\partial x^{2}} \dfrac{\partial ^{2}\Psi _{j} }{ \partial x^{2}}{\rm d}x \end{cases} $ .

2.3.7 Kane equation

For a spacecraft with 18 degrees of freedom, the Kane equation is expressed as:

(49) \begin{align}{F}_{r}^{*}+F_{r}=0\left(r=1,\cdots ,18\right) \end{align}

Substituting Eqs. (35) and (36)into Eq. (49) yields the complete system of dynamic equations for this spacecraft. Rearranging these equations yields the following standard matrix form:

(50) \begin{align}\boldsymbol{M}\, \dot{\boldsymbol{\kern-1pt U}}=\boldsymbol{C}+\boldsymbol{F}+\boldsymbol{D}\boldsymbol{\varLambda } \end{align}

where M is the mass matrix of the spacecraft system’s equations of motion; C is the18 × 1 -dimensional matrix formed by the residual components of the generalized inertial forces after removing the first-order derivatives of the generalized speeds; F is the 18 × 1 -dimensional matrix formed by the residual components of the generalized active forces after removing the Lagrange undetermined multiplier terms; D is the unit vector matrix obtained by normalizing the external normal vectors; $\boldsymbol{\varLambda}$ represents the Lagrangian multiplier matrix, i.e., $\boldsymbol{\varLambda}$ = [λ ₁, λ ₂]^T.

Additionally, the spacecraft body’s attitude relative to the inertial reference frame is described using Euler quaternions, with its kinematic evolution equation given by:

(51) \begin{align} \dot{\boldsymbol{q}} &=\left(1/2\right)S\left(\boldsymbol{q}\right){}^{N}{\boldsymbol{\omega }}{_{}^{B}} \\[-10pt] \nonumber\end{align}

(52) \begin{align} S\left(\boldsymbol{q}\right) &=\boldsymbol{q}_{0}\boldsymbol{E}_{3}+\left[\begin{array}{ll} {\boldsymbol{q}}_{{\rm v}}^{{\rm x}} &\quad \boldsymbol{q}_{{\rm v}}\\[3pt] \boldsymbol{q}_{{\rm v}} &\quad 0 \end{array}\right] \\[10pt] \nonumber\end{align}

where: q = [q ₀, q ₁, q ₂, q ₃]^T represents the attitude four-element vector, q _v = [q ₁, q ₂, q ₃]^T denotes the vector part of the Euler quaternion, and q ₀ represents the scalar part of the Euler quaternion; ^N $\boldsymbol{\omega}$ ^B represents the angular velocity of the spacecraft’s attitude relative to the inertial reference frame in the spacecraft’s principal rigid hub coordinate frame; E ₃ is the unit matrix; ${\boldsymbol{q}_{\textrm{v}}^{\textrm{x}}}$ is the cross product matrix of q _v, defined as:

\begin{align*}{\boldsymbol{q}}_{{\rm v}}^{{\rm x}}=\left[\begin{array}{ccc} 0 &\quad -q_{3} &\quad q_{2}\\[3pt] q_{3} &\quad 0 &\quad -q_{1}\\[3pt] -q_{2} &\quad q_{1} &\quad 0 \end{array}\right] \end{align*}

Therefore, the spacecraft attitude matrix (direction cosine matrix) is represented in Euler quaternions as:

(53) \begin{align}{}^{N}{\boldsymbol{C}}{_{}^{B}}=\left(2{q_{0}}^{2}-1\right)\boldsymbol{E}_{3}+2\boldsymbol{q}_{v}{\boldsymbol{q}_{v}}^{{\rm T}}-2q_{0}{\boldsymbol{q}}_{{\rm v}}^{{\rm x}} \end{align}

Finally, combining Eqs. (44), (50) and (51) yields the coupled system of dynamic equations for the spacecraft:

(54) \begin{align} \begin{cases} \dot{\boldsymbol{q}}=\left(1/2\right)S\left(\boldsymbol{q}\right){}^{N}{\boldsymbol{\omega }}{_{}^{B}}\\[3pt] \boldsymbol{\varLambda }=-\left(\boldsymbol{A}\boldsymbol{M}^{-1}\boldsymbol{D}\right)^{-1}\left[\boldsymbol{A}\boldsymbol{M}^{-1}\left(\boldsymbol{C}+\boldsymbol{F}\right)+\boldsymbol{B}\right]\\[3pt] \, \dot{\boldsymbol{\kern-1pt U}}=\boldsymbol{M}^{-1}\left(\boldsymbol{C}+\boldsymbol{F}+\boldsymbol{D}\boldsymbol{\varLambda }\right) \end{cases} \end{align}