# Dynamic Analysis of Tethered Satellites with a Payload Moving Along a Flexible Tether
## 3 Discretized Equations of Motion
By assembling all the element matrices and vectors, the following global matrix-vector equation can be obtained:  
\[K=\frac{1}{4}\sum_{e=1}^{N}\rho A\dot{\mathbf{d}}_{e}^{\mathrm{T}}\left(\int_ {-1}^{1}\mathbf{N}^{\mathrm{T}}\mathbf{N}\mathrm{d}_{e}^{x}\right)\dot{ \mathbf{d}}_{e}+\frac{1}{2}m_{A}\dot{\mathbf{d}}_{1}^{\mathrm{T}}\left( \mathbf{N}^{\mathrm{T}}\mathbf{N}\right)_{\xi=-1}\dot{\mathbf{d}}_{1}\]3  
where \(\mathbf{M}\) is the global mass matrix, \(\mathbf{G}\) is the global gyroscopic matrix, \(\mathbf{K}\) is the global stiffness matrix, \(\mathbf{q}\) is the global nonlinear force vector, and \(\mathbf{d}\) is the global nodal coordinate vector, given by  
\[K=\frac{1}{4}\sum_{e=1}^{N}\rho A\dot{\mathbf{d}}_{e}^{\mathrm{T}}\left(\int_ {-1}^{1}\mathbf{N}^{\mathrm{T}}\mathbf{N}\mathrm{d}_{e}^{x}\right)\dot{ \mathbf{d}}_{e}+\frac{1}{2}m_{A}\dot{\mathbf{d}}_{1}^{\mathrm{T}}\left( \mathbf{N}^{\mathrm{T}}\mathbf{N}\right)_{\xi=-1}\dot{\mathbf{d}}_{1}\]4  
The force vector \(\mathbf{q}\) is a nonlinear function of the coordinate vector \(\mathbf{d}\). The size of the matrices was \((4N+4)\times(4N+4)\), and the size of the vectors was \((4N+4)\times 1\). To obtain the dynamic responses of the tethered satellite system with a moving payload, the Newmark time integration method [27] was applied to the nonlinear matrix-vector equation in (27). The Newmark method is a second-order accurate, implicit time integration method that permits larger time step sizes compared to explicit methods like the Runge-Kutta method. It is particularly well-suited for long-duration problems, such as the satellite behavior analysis in this study. The Newmark method can experience convergence issues in time-response calculations dominated by the high-frequency modes of flexible bodies, such as in impact analysis. However, it remains robust against numerical instability when applied to satellite response calculations, where high-frequency modes are not predominant. Before applying the Newmark method, the displacement, velocity, and acceleration vectors at time \(t_{n}\), denoted as \(\mathbf{d}(t_{n})\), \(\mathbf{\dot{d}}(t_{n})\), and \(\mathbf{\dot{d}}(t_{n})\), are approximated as \(\mathbf{d}_{n}\), \(\mathbf{v}_{n}\), and \(\mathbf{a}_{n}\), respectively. Denoting \(t_{n+1}=t_{n}+\Delta t\), where \(\Delta t\) is the time step size, the global matrix-vector equation at time \(t_{n+1}\) can be approximated as  
\[K=\frac{1}{4}\sum_{e=1}^{N}\rho A\dot{\mathbf{d}}_{e}^{\mathrm{T}}\left(\int_ {-1}^{1}\mathbf{N}^{\mathrm{T}}\mathbf{N}\mathrm{d}_{e}^{x}\right)\dot{ \mathbf{d}}_{e}+\frac{1}{2}m_{A}\dot{\mathbf{d}}_{1}^{\mathrm{T}}\left( \mathbf{N}^{\mathrm{T}}\mathbf{N}\right)_{\xi=-1}\dot{\mathbf{d}}_{1}\]5  
which is known as the balance equation. In the Newmark time integration method, the displacement and velocity vectors at time \(t_{n+1}\), denoted as \(\mathbf{d}_{n+1}\) and \(\mathbf{v}_{n+1}\), are given by  
\[K=\frac{1}{4}\sum_{e=1}^{N}\rho A\dot{\mathbf{d}}_{e}^{\mathrm{T}}\left(\int_ {-1}^{1}\mathbf{N}^{\mathrm{T}}\mathbf{N}\mathrm{d}_{e}^{x}\right)\dot{ \mathbf{d}}_{e}+\frac{1}{2}m_{A}\dot{\mathbf{d}}_{1}^{\mathrm{T}}\left( \mathbf{N}^{\mathrm{T}}\mathbf{N}\right)_{\xi=-1}\dot{\mathbf{d}}_{1}\]6  
\[K=\frac{1}{4}\sum_{e=1}^{N}\rho A\dot{\mathbf{d}}_{e}^{\mathrm{T}}\left(\int_ {-1}^{1}\mathbf{N}^{\mathrm{T}}\mathbf{N}\mathrm{d}_{e}^{x}\right)\dot{ \mathbf{d}}_{e}+\frac{1}{2}m_{A}\dot{\mathbf{d}}_{1}^{\mathrm{T}}\left( \mathbf{N}^{\mathrm{T}}\mathbf{N}\right)_{\xi=-1}\dot{\mathbf{d}}_{1}\]7  
where \(\beta\) and \(\gamma\) are the algorithmic parameters given by \(1/4\) and \(1/2\), respectively. Using the displacement, velocity, and acceleration vectors at time \(t_{n}\), these vectors were updated at time \(t_{n+1}\).