Add flexible structure model

[200km]

Overview

Implement a flexible structure vibration model.

Details

For detailed attitude control analysis, flexible structures such as solar paddles should be modeled.

Conditions for close

When the flexible structure model is implemented.

Supplement

NA

Note

• Write a conclusion when closing the issue.

[200km]

方針

• S2Eの中で有限要素法を解くなどはせず、外部の解析ツールから出力された構造パラメータをS2Eに入力して、S2E側ではダイナミクスを解くだけ
• 構造パラメータの定義も整理していく
• とりあえず片翼SAPで進めて良いが、最終的には両翼SAPに対応できるようにしたい

今後の流れ

• @t-hosonuma から @fujiy へ、参考になる数式や参考文献の情報を渡す
• @fujiy　の方でトルク形式の数式に変換する
  • 適宜、議論しながら
• 問題なさそうなら実装に進む

[fujiy]

いただいた資料を参考に考えてみたものを一旦まとめます：

参考文献：

• https://stage.tksc.jaxa.jp/taurus/member/miyazaki/old/publication/M_2009_mukumoto_yoshihiro.pdf
• https://stage.tksc.jaxa.jp/taurus/member/miyazaki/old/publication/B_2015_yoshihara_kai.pdf

運動方程式

衛星本体の剛体運動と， $n$ 個の固有モードに分解された柔軟構造物との連成運動を考える． 柔軟構造物の固有モードについては拘束モードモデルを考え，予め解析して算出しておく（衛星全体でモード解析し，本体の並進・回転モードも一緒に扱えば非拘束モードになる？）．

パラメータとして $m \in \mathbb{R}$：衛星全体の質量 $I \in \mathbb{R}^{3\times3}$：衛星全体の慣性テンソル $\delta_0 \in \mathbb{R}^{3 \times n}$：零次干渉行列（本体の並進運動と各モードとの結合係数） $\delta_1 \in \mathbb{R}^{3 \times n}$：一次干渉行列（本体の回転運動と各モードとの結合係数） $M \in \mathbb{R}^{n \times n}$：柔軟モードの質量行列（対角行列） $C \in \mathbb{R}^{n \times n}$：柔軟モードの減衰行列（対角行列） $K \in \mathbb{R}^{n \times n}$：柔軟モードの剛性行列（対角行列） が与えられているものとする．

$x \in \mathbb{R}^3$：本体の変位 $\theta \in \mathbb{R}^3$：本体の回転角 $u \in \mathbb{R}^n$：柔軟モードの変位 $F \in \mathbb{R}^3$：本体に加わる外力 $T \in \mathbb{R}^3$：本体に加わる外力トルク

とすると，運動方程式は

1. 本体の並進運動： $m \ddot{x} + \delta_0 \ddot{u} = F$
2. 本体の回転運動： $I \ddot{\theta} + \delta_1 \ddot{u} = T$
3. 柔軟モードの運動： $M \ddot{u} + C \dot{u} + K u + \delta_0^T \ddot{x} + \delta_1^T \ddot{\theta} = 0$

と表される．

式(1)(2)を(3)に代入すると

$$ M \ddot{u} + C \dot{u} + K u + \frac{1}{m} \delta_0^T (F - \delta_0 \ddot{u}) + \delta_1^T I^{-1} (T - \delta_1 \ddot{u}) = 0 \tag{4} $$

という線型方程式が得られる．

計算の流れ

時間ステップ毎に，

1. 柔軟モード以外による外力・トルクを全て計算する
2. 衛星全体の質量と慣性テンソル，本体に加わる外力とトルクを取得する
3. 式(4)にこれらを代入し，モード変位の加速度 $\ddot{u}$ について解く
4. 3を用いて，柔軟モードの数値積分を行う
5. 3を用いて，本体に柔軟モードによる並進力 $\delta_0 \ddot{u}$ およびトルク $\delta_1 \ddot{u}$ を加算する
6. 本体の数値積分を行う

とすると，ダイナミクスの計算でコンポーネント毎に発生する力を合算しているこれまでの流れをあまり崩さずに，柔軟モードをコンポーネントの一つとして入れられそうです．

[fujiy]

上の式は座標変換の扱いが怪しいので，その点はもう少し考えてみます

[fujiy]

上の式を修正します：

衛星全体（柔軟構造も含む）について， $m \in \mathbb{R}$：質量 [kg] $I \in \mathbb{R}^{3\times3}$：慣性テンソル [kg・m2] $x \in \mathbb{R}^3$：位置 [m] $\theta \in \mathbb{R}^3$：回転角 [rad] $F \in \mathbb{R}^3$：外力 [N] $T \in \mathbb{R}^3$：外力トルク [Nm]

柔軟構造の $i$ 番目のモードについて， $r \in \mathbb{R}^3$：取り付け位置（Body座標） [m] $\beta_i \in \mathbb{R}^3$：並進運動の刺激係数 [√kg] $\gamma_i \in \mathbb{R}^3$：回転運動の刺激係数 [√kg・rad/m] $\omega_i \in \mathbb{R}$：固有角振動数 [rad/s] $\zeta_i \in \mathbb{R}$：減衰比 [-] $q_i \in \mathbb{R}$：モード質量で正規化した変位 [m・√kg]

運動方程式は，

1. 本体の並進運動： $m \ddot{x} + \sum \beta_i \cdot \ddot{q}_i = F$
2. 本体の回転運動： $I \ddot{\theta} + \sum \gamma_i \cdot \ddot{q}_i = T$
3. $i$番目の柔軟モード： $\ddot{q}_i + 2 \zeta_i \omega_i \dot{q}_i + \omega_i^2 q_i + \beta_i \cdot (\ddot{x} + \dot{\theta} \times (\dot{\theta} \times r) + \ddot{\theta} \times r) + \gamma_i \cdot \ddot{\theta} = 0$

式(1)(2)を3に代入して連立方程式を解いて $\ddot{q}_i$ を求め，改めて本体に柔軟モードの運動による力およびトルクを加える．

（式(3)の柔軟モードの運動方程式は回転するBody座標系から見たものなので，遠心力およびオイラー力を考慮するように修正しています．ただし遠心力はモード解析の時点で考慮されないのでここでも不要かもしれません．柔軟構造の変位は微小としているのでコリオリ力は省きました．）

このような方針でいかかでしょうか？ @200km

[200km]

@fujiy ありがとうございます！大きな方針は良いように思いますが、私もこの部分を専門で扱っているわけではないので、実際には実装して簡単なモデルで検証を行いつつ勧めていくことになるかと思います。
研究室の他のメンバーにも問題なさそうか確認してみます。

[TomokiMochizuki]

@.fujiyの代わりに柔軟構造モデルを構築することになりました。

上の方針では $\beta_i, \gamma_i$ などのパラメータを初期値として入力する必要が生じるが、各衛星のこれらの値を正確に取得することは難しいと思っています。

そういう意味では、バネ・マス・ダンパ系のようなシンプルなモデルにしてしまって良い気がします。

例えばですけど、以下の論文の6章に書かれているようなダイナミクスとかはいかがでしょう？ https://arc.aiaa.org/doi/abs/10.2514/6.2011-6592

(I_b + I_p)\ddot{\psi}_b + I_p \ddot{\psi}_s + I_p \ddot{\psi}_d = T_c

I_p(\ddot{\psi}_b + \ddot{\psi}_s + \ddot{\psi}_d) + c_p\dot{\psi}_d + k_p\psi_d = 0 

《変数定義》 $I_b$: 衛星構体（ブームを除く）の慣性モーメント $I_p$: ブームの慣性モーメント $\psi_b$: 衛星の慣性座標における姿勢角 $\psi_s$: ブームの重心の曲げ角 $\psi_d$: ブームの曲げ角 $T_c$: 衛星に作用するトルク $c_p$: ブーム振動の減衰定数 $k_p$: ブーム振動のバネ定数

[TomokiMochizuki]

上式から $\psi_b$ を消去すると、

\frac{I_b I_p}{I_b + I_p}\left(\ddot{\psi}_s + \ddot{\psi}_d \right)
 + \frac{I_p}{I_b + I_p} T_c + c_p \dot{\psi}_d + k_p \psi_d = 0

となる。

[suzuki-toshihir0]

@200km @t-hosonuma

一応補足しておくと、これまで止まっていた柔軟構造の実装を望月くんに進めてもらおうとしています。内容は上に書いてもらっているとおり、岩田さんの論文（ https://arc.aiaa.org/doi/abs/10.2514/6.2011-6592 ）ベースで簡易モデルを実装してもらうつもりで、既にISSL内でも使ったことがある手法だと思うので特に問題はないかな、と思っています。もし何かあればコメントいただければ嬉しいです。