このページは機械翻訳を使用して翻訳されました。最新版の英語を参照するには、ここをクリックします。

Spacecraft Dynamics ブロックを使い始める

この例では、Aerospace Blockset ™ の Spacecraft Dynamics ブロックを使用して、宇宙機または宇宙機群の 6 自由度の剛体ダイナミクスをモデル化する方法を示します。

Spacecraft Dynamics ブロックは、数値積分を使用して宇宙機の並進および回転ダイナミクスをモデル化します。時間の経過に伴う 1 つまたは複数の宇宙機の位置、速度、姿勢、角速度を計算します。最も正確な結果を得るには、許容誤差設定が低い (1e-8 未満) 可変ステップ ソルバーを使用します。ミッションの要件に応じて、より大きな許容値を使用して速度を上げることができます。そうすると、ソリューションの精度に影響する可能性があります。

軌道状態を軌道要素のセットとして、または位置と速度の状態ベクトルとして定義します。軌道状態を伝播するために、ブロックは現在の中心天体に対して選択された重力モデルを使用します。ブロックには、ブロックへの入力として提供される外部加速度と力も含まれます。

姿勢状態は、四元数、方向余弦行列 (DCM)、またはオイラー角を使用して定義されます。姿勢状態を伝播するために、ブロックはブロックへの入力として提供されるモーメントと、ブロック上で定義された質量プロパティを使用します。

このドキュメントでは、ブロックで使用できるさまざまなオプションと構成について説明し、宇宙機のコンスタレーションをモデル化する方法を説明し、低軌道観測衛星用の Spacecraft Dynamics ブロックを実装する Simulink ® モデルの例を示します。最後に、ブロックで使用される方程式が提示されます。

ブロックの説明

Spacecraft Dynamics ブロックは、Simulink ライブラリ ブラウザー (Aerospace Blockset$\to$Spacecraft $\to$Spacecraft Dynamics) にあります。または、Simulink モデル キャンバスのクイック挿入ダイアログに「Spacecraft Dynamics」と入力して見つけることもできます。このセクションでは、Simulink プロパティ インスペクター (モデリング タブの デザイン の下) から表示される、ブロックで使用可能なオプションの概要を示します。

メインタブ

メイン タブには、ブロックレベルの構成パラメータが含まれています。このタブのすべてのパラメータは、ブロックで定義されているすべての宇宙機に適用されます。

以下を含めるかどうかを指定できます。

• 物体力はボディ座標系で定義されます

• ボディモーメントはボディ座標系ポートで定義されます

• 外部加速度は、ブロックの内部計算には含まれない、軌道伝播における摂動加速度を含めます。デフォルトでは、ブロックは中心天体の重力を計算して軌道伝播に使用します (以下の 中心天体 タブ を参照)。伝播に含めることができる追加の摂動加速度の例としては、大気抵抗、第三物体の重力、太陽放射圧によるものがあります。外部加速度座標フレームに設定された値に応じて、慣性 (ICRF) 座標系または固定フレーム座標系で摂動加速度を提供できます。各中央ボディに使用される固定フレーム座標系の詳細については、ブロック参照ページの「座標系」セクションを参照してください。

状態ベクトル出力座標フレームは、ブロックからの位置と速度の状態出力が慣性 (ICRF) 座標系か固定フレーム座標系かを制御します。

常に慣性フレーム内にあるブロックから合計慣性加速度を出力するかどうかも指定できます。この値は、内部で計算された中心物体の重力と、ブロックに入力として提供される物体力と外部加速度の寄与を含む合計加速度です。注意: 加速出力ポートは診断目的のみに使用されます。この信号を入力としてブロックにフィードバックすることは有効なワークフローではありません。

開始日時は、Simulink モデルの開始時刻 ${\mathit{t}}_0$ に対応する最初の日時です。これは、ブロックで提供されるすべての初期条件の想定されるエポックです。オプションで、現在の日付/時刻を出力を選択して、ブロックから時間信号を出力し、シミュレーションの他の場所で使用することができます。

質量タブ

宇宙機の質量特性をモデル化するために、Fixed、 Simple variable、Custom Variable の 3 つの質量タイプが利用できます。

質量タイプがFixedの場合、質量と慣性テンソルはシミュレーション全体を通じて質量と慣性テンソルに指定された値で一定に保たれます。質量流量と慣性の変化率はゼロになります。

質量タイプがSimple variableの場合、シミュレーション中に宇宙機の質量特性を変化させるために単純なアプローチが採用されます。

初期の 質量、空の質量 (乾燥)、および 全質量 (湿潤) が定義されます。質量流量は入力ポート (dm/dt) を介してブロックに提供されます。この値は積分され、シミュレーションの各タイムステップでの現在の質量を計算します。

同様に、この構成では、空の宇宙機構成と完全な宇宙機構成の慣性テンソル値を指定します。現在のテンソルは、現在の質量値に基づいて、空の慣性テンソルと完全な慣性テンソルの間の線形補間によって近似されます。

オプションで、ブロックに別の入力ポート (${\mathit{V}}_{{\mathrm{re}}_{\mathit{b}}}$) を追加し、パラメータ 質量流量相対速度を含める を使用して質量流量相対速度を提供することができます。この相対速度はボディ座標系で提供されます。これは、宇宙機から除去される質量や宇宙機に追加される質量による力の寄与を計算するために使用されます。

現在の質量が空の質量値を下回っているか、または満杯の質量値を超えている場合に、ブロックに提供される質量流量を制限するには、パラメーター 質量が空または満杯の場合に質量流量を制限する を使用します。

最後に、現在の質量に基づいて、現在の燃料の状態をブロック (Fuel Status) から出力できます。現在の質量が 全質量 を超える場合、報告されるステータスは 1 になります。現在の質量が 空の質量 を下回る場合、ステータスは -1 になります。質量が指定された動作範囲内にある場合、ステータスは 0 になります。

質量タイプがCustom variableの場合、宇宙機の質量特性が時間の経過とともにどのように変化するかに関して、より柔軟性が提供されます。ただし、これを行うには、ブロックの外部でさらに多くの値を計算する必要があります。

この構成では、現在の質量 (m)、現在の慣性テンソル (I)、および慣性テンソルの現在の変化率 (dI/dt) のブロック入力ポートが追加されます。

質量流量の相対速度を提供するには、オプションで、質量流量の相対速度を含めるパラメータを使用して、ブロックに別の入力ポート (${\mathit{V}}_{{\mathrm{re}}_{\mathit{b}}}$) を追加できます。この相対速度はボディ座標系で提供されます。これは、宇宙機から除去される質量や宇宙機に追加される質量による力の寄与を計算するために使用されます。したがって、このパラメータを有効にすると、ブロックに質量流量を提供するための追加ポートが追加されます (dm/dt)。

軌道タブ

軌道タブでは、初期状態形式(Orbital elements、ICRF state vector、またはFixed-frame state vector)の値に応じて、宇宙機の初期条件を軌道要素のセットとして、または位置と速度の状態ベクトルとして定義します。

状態ベクトル オプションの場合、指定された座標フレームの 開始日時 に対応する初期位置と速度を指定します。

初期状態形式 オプション Orbital elements は、パラメータ 軌道タイプ によってさらに分解されます。

軌道タイプが Keplerian の場合、6 つのケプラー軌道要素の従来のセットを指定します。

• 半長径 ($\mathit{a}$)

• 偏心 ($\mathit{e}$)

• 傾向 ($\mathit{i}$)

• 昇交点の赤経 - RAAN ($\Omega$)

• 近点引数 ($\omega$)

• 真の異常 ($\nu$)。

軌道要素を指定する場合、次の 3 つの軌道タイプによって未定義の要素が生成されます。

• 軌道が赤道軌道（傾斜角がゼロ）の場合、RAAN は定義されません。

• 軌道が円形（離心率がゼロ）の場合、近点引数と真近点角は定義されません。

• 軌道が円形で、かつ赤道軌道の場合、3 つの要素はすべて未定義になります。

これらの条件のモデル化を支援するために、ブロックには Keplerian に加えて、軌道タイプ の 3 つの追加オプションが用意されています。

非円形 (楕円形) の赤道軌道の場合、傾斜角は常にゼロになり、RAAN と 近点引数 は 近点経度 に置き換えられます。近点経度は、ICRF X 軸 (I) と近点の間の角度です。これは、RAAN ($\Omega$) と近点引数 ($\omega$) の合計に等しくなります。

円形の傾斜軌道（非赤道軌道）の場合、離心率は常にゼロになり、近点引数と真近点角は緯度引数に置き換えられます。緯度引数は昇交点と衛星位置ベクトル間の角度です。これは、真近点角 ($\nu$) と近点引数 ($\omega$) の合計に等しくなります。

最後に、円赤道軌道の場合、傾斜角と離心率は常にゼロになり、RAAN、近点引数、真近点角は真経度に置き換えられます。真経度は、ICRF X 軸 (I) と宇宙機の位置ベクトルの間の角度です。これは、真日点角 ($\nu$)、近点引数 ($\omega$)、および RAAN ($\Omega$) の合計に等しくなります。

態度タブ

姿勢タブでは、モデル化される宇宙機の姿勢の初期条件を定義します。パラメータ 姿勢参照座標フレーム を使用すると、慣性 (ICRF) フレーム、固定フレーム、北東下 (NED) フレーム、またはローカル垂直ローカル水平 (LVLH) フレームに対する姿勢を定義できます。ブロックに提供される初期姿勢と本体の角速度パラメータは、指定されたフレームに対して定義されているものと想定されます。ブロックからの姿勢とボディの角速度の出力もこのフレームを使用します。

姿勢に使用する表現方法を指定するには、姿勢表現 パラメータを使用します。選択した値に応じて、初期姿勢 パラメータは 初期クォータニオン、初期 DCM、または 初期オイラー角 として表示されます。提供された初期条件の次元がその表現と一致することが期待されます。さらに、ブロックからの姿勢出力ポートは指定された表現を使用します。

初期姿勢変化率、初期ボディ角速度 PQR、および対応する出力ポート ($\omega$) は、姿勢表現 の選択に関係なく、常に角速度として定義されます。

ブロックから 総慣性角加速度 ($\dot{\omega }$) を出力するかどうかも指定できます。この出力は常に慣性 (ICRF) フレームに対して定義されます。重力勾配トルクを含める が選択されている場合、この値はブロックへの入力として提供されるモーメントと内部で計算された重力勾配トルクによる合計角加速度です。注意: 角加速度出力ポートは診断目的のみに使用されます。この信号を入力としてブロックにフィードバックすることは有効なワークフローではありません。

有効にすると、重力勾配トルクの計算では中心物体が球体として扱われます。重力勾配トルクによる全体的な寄与は小さいです。一般に、中心体を球体として扱うことは、ほとんどのアプリケーションでは十分です。より高い精度が必要な場合は、重力勾配トルク値を外部で計算し、モーメントとしてブロックに提供することができます。ブロックによって実装される方程式については、以下のブロック方程式のセクションを参照してください。

中央本文タブ

中心天体タブを使用して、宇宙機が周回する天体の物理的特性、重力ポテンシャル モデル、および方向に関する情報を提供します。地球、月（ルナ）だけでなく、太陽系のすべての惑星も利用可能です。カスタム中央ボディも定義できます。まず、パラメータ 重力ポテンシャル モデル のさまざまなオプション (None、Point-mass、Oblate ellipsoid (J2)、および Spherical harmonics) を確認します。

None はすべての中心体で利用可能です。このオプションでは、システム方程式に内部的に計算された重力加速度は含まれません。独自の重力モデルを使用する場合は、このオプションを外部加速度入力ポートと組み合わせて使用します。オプション None をカスタム中心本体とともに使用する場合は、惑星の 回転速度 のみが必要です。

Point-mass はすべての中心体で利用可能です。このオプションは、中心物体を質点として扱い、ニュートンの万有引力の法則を使用して重力加速度を計算します。オプション Point-mass をカスタム中心天体とともに使用する場合は、赤道半径、平坦化、重力パラメーター ($\mu$)、および惑星の 回転速度 を指定する必要があります。

Oblate ellipsoid (J2) はすべての中心体で利用可能です。このオプションには、中心天体の扁平度を考慮した、2 次帯状調和重力係数 ${\mathit{J}}_2$ の摂動効果が含まれます。オプション Oblate ellipsoid (J2) をカスタム中心天体とともに使用する場合は、赤道半径、平坦化、重力パラメーター ($\mu$)、第 2 次ゾーン調和関数 (J2)、および惑星の 自転速度 を指定する必要があります。

Spherical harmonics は、中央本体が Earth、Moon、Mars、または Custom に設定されている場合にのみ使用できます。各中心体で使用可能な 球面調和モデル オプションが次の表にリストされています。

また、選択した球面調和関数モデルでサポートされている最大度数より低い値を Degree に指定する必要があります。各モデルの推奨度数と最大度数の値は以下の通りです。

オプション Spherical harmonics をカスタム中心体とともに使用する場合は、惑星の 回転速度、球面調和係数ファイル (.mat)、および 度 を指定する必要があります。このファイルの詳細については、Spacecraft Dynamics ブロックのリファレンス ページのパラメータの説明を参照してください。

ブロックで使用されるすべての惑星定数は、NASA JPL Planetary and Lunar Ephemerides DE405 からのものです。

すべての J2 定数値は、NASA 宇宙科学データコーディネートアーカイブ (NSSDCA) から取得されます。

代替の定数値が必要な場合は、Central body. に Custom オプションを使用します。

重力に加えて、中心物体 タブには中心物体の方向に関する情報が含まれます。使用可能なパラメータは、現在選択されている Central ボディによって異なります。 Earth、Moon、Custom を除くすべての中心天体は、 地図座標と回転要素に関する IAU/IAG ワーキング グループの報告書からの惑星回転極と子午線の定義を使用します。2006 .地球、月、カスタムに固有のオプションについては以下で説明します。

中心体がEarthの場合、ブロックが使用する固定フレーム座標系はITRFです。デフォルトでは、ICRF と ITRF 間の変換では、パラメータ IERS EOP データ ファイル に提供される地球方向パラメータ (EOP) データが使用されます。最新の EOP データ ファイルを生成するには、Aerospace Toolbox 関数 aeroReadIERSData() を使用します。この関数は、IERS データ サーバーを呼び出し、最新の EOP データを MAT ファイルに保存します。変換から地球方向パラメータ データを除外するには、地球方向パラメータ (EOP) を使用する をオフにします。

中心天体がMoonの場合、月の秤動角を入力を選択して、月の秤動角をブロックへの入力として提供できます。このオプションを選択すると、ブロックに入力ポートが追加されます。この場合、シミュレーションの各タイム ステップで秤動角が提供され、固定フレームと ICRF 間の変換に使用されます。Moon Libration ブロックを使用して、秤動角を計算できます。秤動角を使用する場合、Moon の固定フレーム座標系は平均地球/極軸フレーム (ME) です。このフレームは 2 つの変換によって実現されます。まず、ブロックは ICRF フレームの値を主軸システム (PA) に変換します。主軸システムは、ブロックへの入力として提供される秤動角によって定義される軸です。次に、ブロックは、地図座標と回転要素に関する IAU/IAG ワーキング グループのレポート からの固定回転を使用して、状態を ME システムに変換します。2006 .月の秤動角を入力 チェックボックスをオフにすると、固定フレームは、地図座標と回転要素に関する IAU/IAG ワーキンググループの レポートで定義されている回転極と本初子午線の方向によって定義されます。2006 .

中心体がCustomの場合、パラメータ中心体スピン軸ソースの値に応じて、ブロックに回転極と子午線データを提供するためのオプションが 2 つあります。各タイムステップで現在の赤経、赤緯、および回転速度の値をブロックへの入力として提供するには、ソースを Port に設定します。J2000 (JD 2451545.0、つまり 2000 年 1 月 1 日 12 時間 TDB) での赤経、赤緯、および回転角の初期条件と、各値に対応する変化率を提供するには、ソースを Dialog に設定します。これらのパラメータは、地図座標と回転要素に関する IAU/IAG ワーキング グループの レポートに記載されている用語と方程式と一致しています。2006 .

最後に、すべての中心体に対して、出力クォータニオン (ICRF から固定フレーム) を選択して、ICRF から固定フレームへの位置変換を実行するクォータニオンをオプションで出力できます。

単位タブ

単位 タブでは、ブロックで使用される単位系、角度単位 (Degrees または Radians)、および時間形式を定義します。

時間形式がJulian dateの場合、開始日時とブロックのオプションの時間出力ポートはスカラーのユリウス日付値を使用します。Gregorian に設定すると、両方の値は [年、月、日、時間、分、秒] 形式の 1x6 配列になります。パラメータ Units の各オプションに対応する単位を以下の表に示します。単位 が変更されると、ブロック上の各パラメータとポート ラベルの予想される単位が自動的に更新されます。

衛星コンスタレーションのモデリング

ここまで、Spacecraft Dynamics ブロックを使用して単一の宇宙機をモデル化しました。ただし、ブロックは衛星/宇宙機のコンスタレーションをモデル化するように構成することもできます。モデル化される宇宙機の数は、提供される初期条件のサイズによって決まります。質量、軌道、または 姿勢 タブのパラメータに複数の値が指定されている場合、ブロックは衛星のコンスタレーションを出力します。単一の値が指定されたパラメータは拡張され、コンスタレーション内のすべての衛星に適用されます。たとえば、6 つの値を含む True anomaly を除くブロックのすべてのパラメータに 1 つの値が指定されている場合、True anomaly のみが変化する 6 つの衛星のコンスタレーションが作成されます。

この動作は、すべての宇宙機の初期条件（質量、軌道、または 姿勢 タブのすべてのダイアログ ボックス）に適用されます。上記の 質量、軌道、または 姿勢 タブの初期条件パラメータには、すべての衛星に拡張された単一の値、または衛星ごとに 1 つずつ、6 つの個別の値が含まれている必要があります。

同じ拡張動作はブロック入力ポートにも適用されます。すべての入力ポートは、拡張を 期待 します 月の秤動角 $\phi \theta \psi$ (中心天体 が Moon の場合) とスピン軸の 赤経、赤緯、および回転角 $\alpha \delta \mathit{W}$ (中心天体 が Custom の場合) をサポートします。月の秤動角と自転軸の方向の入力は時間に依存する値であるため、モデル化されているすべての宇宙機に常に適用されます。他のすべてのポートは、モデル化されているすべての宇宙機に拡張された単一の値、または各宇宙機に適用される個別の値 (上記の例では 6) を受け入れます。

月の軌道のモデル化

このポート拡張動作を示すために、軌道に沿って 200 km 離れた双子の月周回衛星が存在する新しいシナリオを考えてみましょう。各衛星は独立して動作するため、それぞれに異なる力とモーメントが適用されます。しかし、地球の重力の影響を両方の衛星に対する摂動加速度として考慮に入れたいと考えています。200km にわたる月の軌道における地球の重力加速度の差は無視できると仮定します。結果として得られたブロックを以下に示します。

各衛星には個別の力とモーメントの入力値がありますが、単一の外部加速度入力が拡張され、両方の宇宙機に適用されます。前述のように、月の秤動角 $\phi \theta \psi$ は常に宇宙機に依存しません。

ブロックからの状態出力は、モデル化されている宇宙機の合計数と常に一致し、行は個々の宇宙機に対応します。ブロックからは、現在の時間 ${\mathit{t}}_{\mathrm{utc}}$ と慣性フレームから固定フレームへの変換 ${\mathit{q}}_{\mathrm{icrf2ff}}$ という 2 つの時間依存出力もあります。

Simulink モデル例

ここで、Spacecraft Dynamics ブロックを使用して地球観測衛星をシミュレートするサンプルモデルを調べてみましょう。

mdl = "SpacecraftDynamicsBlockExampleModel";
open_system(mdl);

衛星は高度約500kmのほぼ円形の低地球軌道（LEO）にあります。軌道伝播には、次数を 120 に設定した地球球面調和関数モデル EGM2008 を使用します。Aerospace Toolbox に含まれるデフォルト ファイル aeroiersdata.mat の地球方向パラメータ データを使用します。衛星の質量特性は固定されており、質量は 1kg 、単純慣性テンソルは $\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]$ です。ミッション開始日は2020年1月1日12:00:00で、実行時間は6時間です。

姿勢制御を行うために、Aerospace Blockset の Attitude Profile ブロックを、シンプルな PD コントローラーに接続して使用します。望ましい姿勢は、衛星本体のz軸を高度10mの地理座標$\left[{42}^{\circ },{-71}^{\circ }\right]$に合わせることです。2 次制約として、body-x 軸を、ローカル垂直、ローカル水平 (LVLH) フレームの y 軸に揃えます。(ほぼ)円軌道では、LVLH フレームの y 軸は衛星の移動方向を指します。この配置により、衛星は通過するたびに地理的に興味のある地点の上空を移動しながら「前方」を向き続けます。また、Attitude Profile ブロックを使用して、衛星を地球の天底、別の地理的位置、JPL Ephemerides DE405 の天体、またはブロックへの入力として提供される任意のカスタム ベクトルに合わせることもできます。

モデルには、1U CubeSat を視覚化するように構成された Simulink 3D Animation ™ ワールドが含まれています。このブロックは Simulink 3D Animation を必要とするため、デフォルトでコメント アウトされています。視覚化を有効にし、視覚化のプロパティを変更するには、ブロックをダブルクリックします。

このモデルには、スコープに接続された Spacecraft Dynamics ブロックからの線形加速度出力と角加速度出力もあります。これらの出力をシミュレーション ループの一部として使用しないでください。慣性線形加速はシミュレーション全体を通じてスムーズであり、これは並進 (デルタ V) 操作を実行していないため予想されることです。角加速度プロットでは、地理的に興味のある地点を通過するときに、より大きな加速度値と一致することがわかります。関心のあるポイントを直接通過する場合 (最初のパス)、必要な角速度の変化は、関心のあるポイントをより浅い角度で通過する場合 (後続のパス) よりもはるかに大きくなります。

ブロック方程式

ここで、ブロックによって実装される方程式を調べて、ブロックが各タイムステップで出力値を計算する方法をよりよく理解します。

並進システム方程式

移動動作は以下によって制御されます。

$${\overrightarrow{\mathit{a}}}_{\mathrm{icrf}}={\overrightarrow{\mathit{a}}}_{\mathrm{central}\mathrm{body}\mathrm{gravity}}+\mathrm{body2inertial}\left(\frac{\overrightarrow{{\mathit{F}}_{\mathit{b}}}}{\mathit{m}}\right)+{\overrightarrow{\mathit{a}}}_{\mathrm{applied}}$$

$${\overrightarrow{\mathit{a}}}_{\mathrm{icrf}}\underset{\mathrm{integrate}}{⟹}{\overrightarrow{\mathit{v}}}_{\mathrm{icrf}},{\overrightarrow{\mathit{r}}}_{\mathrm{icrf}}$$

ここで、

${\overrightarrow{\mathit{a}}}_{\mathrm{central}\mathrm{body}\mathrm{gravity}}$ は、現在のブロック パラメータの選択に基づく中心体の重力です。

${\overrightarrow{\mathit{a}}}_{\mathrm{applied}}$ は、ブロックの外部加速度入力ポートに提供されるユーザー定義の加速度です。

$\overrightarrow{{\mathit{F}}_{\mathit{b}}}$ は、ICRF フレーム (慣性) に対するボディ座標系の物体力です。

$\mathit{m}$ は現在の宇宙機の質量です (以下の 質量方程式 を参照)。

$\mathrm{body2inertial}\left(\right)$ は回転体固定座標系から慣性 ICRF 座標系への変換であり、その結果、力による加速度の寄与は次のようになります:$\mathrm{body2inertial}\left(\frac{\overrightarrow{{\mathit{F}}_{\mathit{b}}}}{\mathit{m}}\right)={\overrightarrow{\mathit{a}}}_{{\mathrm{icrf}}_{\mathrm{forces}}}=\mathrm{quatrotate}\left({\mathit{q}}_{\mathrm{b2icrf}},{\overrightarrow{\mathit{a}}}_{\mathit{b}}\right)$

ここで、

$${\overrightarrow{\mathit{a}}}_{\mathit{b}}=\frac{{\overrightarrow{\mathit{F}}}_{\mathit{b}}}{\mathit{m}}=\frac{{\overrightarrow{\mathit{F}}}_{{\mathit{b}}_{\mathrm{input}}}+\dot{\mathit{m}}{\overrightarrow{\mathit{v}}}_{\mathrm{re}}}{\mathit{m}}$$

${\overrightarrow{\mathit{r}}}_{\mathit{b}}=\mathrm{quatrotate}\left({\mathit{q}}_{\mathrm{icrf2b}},{\overrightarrow{\mathit{r}}}_{\mathrm{icrf}}\right)$

${\overrightarrow{\mathit{F}}}_{{\mathit{b}}_{\mathrm{input}}}$ は、ブロック本体の力入力ポートに提供される力です。

${\overrightarrow{\mathit{v}}}_{\mathrm{re}}$ は、ボディ座標系においてボディ座標系を基準として質量流 ($\dot{\mathit{m}}$) がボディから排出されるかボディに追加される相対速度です。

${\mathit{q}}_{\mathrm{icrf2b}}$ は、慣性 ICRF フレームに対する物体の受動的な四元数回転です。

${\overrightarrow{\omega }}_{\mathrm{icrf2b}}$ は、慣性 ICRF フレームに対する物体の角速度です。

回転システム方程式

回転運動は以下によって制御されます。

$${\dot{\overrightarrow{\omega }}}_{\mathrm{icrf2b}}=\left[\overrightarrow{{\mathrm{M}}_{\mathit{b}}}-{\overrightarrow{\omega }}_{\mathrm{icrf2b}}\times \left({\mathit{I}}_{\mathit{mom}}{\overrightarrow{\omega }}_{\mathrm{icrf2b}}\right)-{\dot{\mathit{I}}}_{\mathit{mom}}{\overrightarrow{\omega }}_{\mathrm{icrf2b}}\right]\mathrm{inv}\left({\mathit{I}}_{\mathit{mom}}\right)$$

$${\dot{\overrightarrow{\omega }}}_{\mathrm{icrf2b}}\underset{\mathrm{integrate}}{⟹}{\overrightarrow{\omega }}_{\mathrm{icrf2b}},\overrightarrow{{\mathit{q}}_{\mathrm{icrf2b}}}$$

ここで、

${\mathit{I}}_{\mathit{mom}}=\left[\begin{array}{ccc}{\mathit{I}}_{\mathrm{xx}}& -{\mathit{I}}_{\mathrm{xy}}& -{\mathit{I}}_{\mathrm{xz}}\\ -{\mathit{I}}_{\mathrm{yx}}& {\mathit{I}}_{\mathrm{yy}}& -{\mathit{I}}_{\mathrm{yz}}\\ -{\mathit{I}}_{\mathrm{zx}}& -{\mathit{I}}_{\mathrm{zy}}& {\mathit{I}}_{\mathrm{zz}}\end{array}\right]$ は、物体の原点に対する慣性テンソルです (以下の 質量セクション を参照)。

${\dot{\mathit{I}}}_{\mathit{mom}}$ は慣性テンソルの変化率です (以下の 質量セクション を参照)。

$\mathrm{inv}\left(\right)$ は 3x3 行列の逆行列です。

$\overrightarrow{{\mathrm{M}}_{\mathit{b}}}={\overrightarrow{\mathit{M}}}_{{\mathit{b}}_{\mathrm{input}}}+{\mathit{M}}_{{\mathit{b}}_{\mathrm{gravity}\mathrm{gradient}}}$ は、ブロック ボディ モーメント入力ポートに提供された値と内部で計算された重力勾配トルクで構成される合計ボディ モーメントです。

$${\mathit{M}}_{{\mathit{b}}_{\mathrm{gravity}\mathrm{gradient}}}=\frac{3\mu }{{\mathit{r}}_{\mathit{b}}^5}{\overrightarrow{\mathit{r}}}_{\mathit{b}}\times {\mathit{I}}_{\mathit{mom}}{\overrightarrow{\mathit{r}}}_{\mathit{b}}$$

$\mu$ は中心天体の標準的な重力パラメータです。

四元数ベクトルの変化率の積分は次のように計算されます。

$\left[\begin{array}{c}\dot{{\mathit{q}}_0}\\ \dot{{\mathit{q}}_1}\\ \dot{{\mathit{q}}_2}\\ \dot{{\mathit{q}}_3}\end{array}\right]=\left[\begin{array}{cccc}0& {\omega }_{\mathit{b}}\left(1\right)& {\omega }_{\mathit{b}}\left(2\right)& {\omega }_{\mathit{b}}\left(3\right)\\ {-\omega }_{\mathit{b}}\left(1\right)& 0& {-\omega }_{\mathit{b}}\left(3\right)& {\omega }_{\mathit{b}}\left(2\right)\\ {-\omega }_{\mathit{b}}\left(2\right)& {\omega }_{\mathit{b}}\left(3\right)& 0& {-\omega }_{\mathit{b}}\left(1\right)\\ {-\omega }_{\mathit{b}}\left(3\right)& -{\omega }_{\mathit{b}}\left(2\right)& {\omega }_{\mathit{b}}\left(1\right)& 0\end{array}\right]\left[\begin{array}{c}{\mathit{q}}_0\\ {\mathit{q}}_1\\ {\mathit{q}}_2\\ {\mathit{q}}_3\end{array}\right]$

Aerospace Toolbox と Aerospace Blockset は、スカラー優先規則を使用して定義された四元数を使用します。

質量

上記のシステム方程式で使用される質量 $\mathit{m}$、質量流量 $\dot{\mathit{m}}$、慣性テンソル ${\mathit{I}}_{\mathrm{mom}}$、および慣性テンソルの変化率 ${\dot{\mathit{I}}}_{\mathrm{mom}}$ は、質量 タブの現在のパラメータ選択に応じて決定されます。

Fixed

このオプションは、宇宙機を固定質量の剛体としてモデル化します。

$\mathit{m}$ は、質量タブのパラメータ 質量に指定された質量です。

$\dot{\mathit{m}}$ はゼロになります。

${\mathit{I}}_{\mathrm{mom}}$ は、質量タブのパラメータ 慣性テンソルに提供される慣性テンソルです。

${\dot{\mathit{I}}}_{\mathrm{mom}}$ はゼロになります。

Simple variable

このオプションは、宇宙機を単純な可変質量の剛体としてモデル化します。

$\mathit{m}$ は ${\mathit{m}}_{\mathrm{full}}$ と ${\mathit{m}}_{\mathrm{empty}}$ の間に制限される質量であり、 $\dot{\mathit{m}}$ を積分することによって計算されます。

$\dot{\mathit{m}}$ はブロックの dm/dt 入力ポートに提供されます。

$${\mathit{I}}_{\mathrm{mom}}=\frac{{\mathit{I}}_{\mathrm{full}}-{\mathit{I}}_{\mathrm{empty}}}{{\mathit{m}}_{\mathrm{full}}-{\mathit{m}}_{\mathrm{empty}}}\left(\mathit{m}-{\mathit{m}}_{\mathrm{empty}}\right)+{\mathit{I}}_{\mathrm{empty}}$$

$${\dot{\mathit{I}}}_{\mathrm{mom}}=\frac{{\mathit{I}}_{\mathrm{full}}-{\mathit{I}}_{\mathrm{empty}}}{{\mathit{m}}_{\mathrm{full}}-{\mathit{m}}_{\mathrm{empty}}}\dot{\mathit{m}}$$

Custom variable

このオプションは、宇宙機を可変質量の剛体としてモデル化し、最高レベルの構成可能性を提供します。

$\mathit{m}$ はブロックの m 入力ポートに提供されます。

質量流量相対速度を含めるが有効になっている場合、$\dot{\mathit{m}}$はブロックdm/dt入力ポートに提供されます。それ以外の場合、システム方程式では値は必要ありません。

${\mathit{I}}_{\mathrm{mom}}$ はブロックの I 入力ポートに提供されます。

${\dot{\mathit{I}}}_{\mathrm{mom}}$ はブロックの dI/dt 入力ポートに提供されます。

中心体重力

中心物体の重力による加速度 ${\overrightarrow{\mathit{a}}}_{\mathrm{central}\mathrm{body}\mathrm{gravity}}$ は、中心物体 タブの現在のパラメータ選択に応じて計算されます。非球面加速度項 (Oblate ellipsoid (J2) および Spherical harmonics) を含む重力モデルの場合、非球面重力は固定フレーム座標系 (地球の場合は ITRF) で計算されます。ただし、数値積分は常に慣性 ICRF 座標系で実行されます。したがって、各タイムステップで、位置と速度の状態が固定フレームに変換され、非球面重力が固定フレームで計算され、結果として得られる加速度が慣性フレームに変換されます。慣性系では、結果として得られる加速度は他の加速度項と合計され、二重積分されて速度と位置が求められます。

Point-mass (すべての中央機関で利用可能)

このオプションは、ニュートンの万有引力の法則を使用して、球状の重力の影響のみを含め、中心物体を質点として扱います。

$${\overrightarrow{\mathit{a}}}_{\mathrm{central}\mathrm{body}\mathrm{gravity}}=-\frac{\mu }{{{\mathit{r}}_{\mathrm{icrf}}}^2}\frac{{\overrightarrow{\mathit{r}}}_{\mathrm{icrf}}}{{\mathit{r}}_{\mathrm{icrf}}}$$

Oblate ellipsoid (J2) (すべての中央機関で利用可能)

このオプションには、球状の重力に加えて、中心天体の扁平度を説明する第 2 次帯状調和重力係数 ${\mathit{J}}_2$ の摂動効果が含まれます。${\mathit{J}}_2$ は、中心天体の完全な球体からの重力の逸脱の大部分を説明します。

$${\overrightarrow{\mathit{a}}}_{\mathrm{central}\mathrm{body}\mathrm{gravity}}=-\frac{\mu }{{{\mathit{r}}_{\mathrm{icrf}}}^2}\frac{{\overrightarrow{\mathit{r}}}_{\mathrm{icrf}}}{{\mathit{r}}_{\mathrm{icrf}}}+\mathrm{fixed2inertial}\left({\overrightarrow{\mathit{a}}}_{\mathrm{nonspherical}}\right)$$

ここで、

$$\begin{array}{l}{\overrightarrow{\mathit{a}}}_{\mathrm{nonspherical}}=\left\{\left[\frac{1}{\mathit{r}}\frac{\partial }{\partial \mathit{r}}\mathit{U}-\frac{{{\mathit{r}}_{\mathrm{ff}}}_{\mathit{k}}}{{\mathit{r}}^2\sqrt{{{{\mathit{r}}_{\mathrm{ff}}}_{\mathit{i}}}^2+{{{\mathit{r}}_{\mathrm{ff}}}_{\mathit{j}}}^2}}\frac{\partial }{\partial \varphi }\mathit{U}\right]{{\mathit{r}}_{\mathrm{ff}}}_{\mathit{i}}\right\}\mathit{i}\\ +\left\{\left[\frac{1}{\mathit{r}}\frac{\partial }{\partial \mathit{r}}\mathit{U}+\frac{{{\mathit{r}}_{\mathrm{ff}}}_{\mathit{k}}}{{\mathit{r}}^2\sqrt{{{{\mathit{r}}_{\mathrm{ff}}}_{\mathit{i}}}^2+{{{\mathit{r}}_{\mathrm{ff}}}_{\mathit{j}}}^2}}\frac{\partial }{\partial \varphi }\mathit{U}\right]{{\mathit{r}}_{\mathrm{ff}}}_{\mathit{j}}\right\}\mathit{j}\\ +\left\{\frac{1}{\mathit{r}}\left(\frac{\partial }{\partial \mathit{r}}\mathit{U}\right){\mathit{r}}_{\mathit{k}}+\frac{\sqrt{{{{\mathit{r}}_{\mathrm{ff}}}_{\mathit{i}}}^2+{{{\mathit{r}}_{\mathrm{ff}}}_{\mathit{j}}}^2}}{{\mathit{r}}^2}\frac{\partial }{\partial \varphi }\mathit{U}\right\}\mathit{k}\end{array}$$

球座標における偏微分が与えられると：

$\begin{array}{l}\frac{\partial }{\partial \mathit{r}}\mathit{U}=\frac{3\mu }{{\mathit{r}}^2}{\left(\frac{{\mathit{R}}_{\mathrm{cb}}}{\mathit{r}}\right)}^2{\mathit{P}}_{2,0}\left[\sin \left(\varphi \right)\right]{\mathit{J}}_2\\ \frac{\partial }{\partial \varphi }\mathit{U}=-\frac{\mu }{\mathit{r}}{\left(\frac{{\mathit{R}}_{\mathrm{cb}}}{\mathit{r}}\right)}^2{\mathit{P}}_{2,1}\left[\sin \left(\varphi \right)\right]{\mathit{J}}_2\end{array}$

ここで、

$\varphi$と$\lambda$は衛星の地心緯度と経度です

${\mathit{P}}_{2,0}$と${\mathit{P}}_{2,1}$はルジャンドル関数である。

${\mathit{R}}_{\mathrm{cb}}$は中心天体の赤道半径である

$\mathrm{fixed2inertial}\left(\right)$ は、遠心加速度とコリオリ加速度を考慮して、固定フレームの位置、速度、加速度を中心物体の中心を原点とする ICRF 座標系に変換します。各中心ボディに使用される固定座標系と慣性座標系の詳細については、Spacecraft Dynamics ブロック参照ページの「座標フレーム」セクションを参照してください。地球に使用される固定フレーム座標フレームは ITRF です。

Spherical Harmonics (地球、月、火星、カスタムで利用可能)

このオプションは、ゾーン、セクター、およびテッセラル高調波を考慮した高次の摂動効果を含めることで、忠実度を高めます。参考までに、2 次ゼロ次ゾーン高調波 ${\mathit{J}}_2$ は $-{\mathit{C}}_{2,0}$ です。球面調和関数モデルは、中心体とジオポテンシャル モデルによって変化する最大次数 $\mathit{l}={\mathit{l}}_{\max }$ までの調和関数を考慮します。

$${\overrightarrow{\mathit{a}}}_{\mathrm{central}\mathrm{body}\mathrm{gravity}}=-\frac{\mu }{{\mathit{r}}^2}\frac{\overrightarrow{{\mathit{r}}_{\mathrm{icrf}}}}{\mathit{r}}+\mathrm{fixed2inertial}\left({\overrightarrow{\mathit{a}}}_{\mathrm{nonspherical}}\right)$$

ここで、

$$\begin{array}{l}{\overrightarrow{\mathit{a}}}_{\mathrm{nonspherical}}=\left\{\left[\frac{1}{\mathit{r}}\frac{\partial }{\partial \mathit{r}}\mathit{U}-\frac{{{\mathit{r}}_{\mathrm{ff}}}_{\mathit{k}}}{{\mathit{r}}^2\sqrt{{{{\mathit{r}}_{\mathrm{ff}}}_{\mathit{i}}}^2+{{{\mathit{r}}_{\mathrm{ff}}}_{\mathit{j}}}^2}}\frac{\partial }{\partial \varphi }\mathit{U}\right]{{\mathit{r}}_{\mathrm{ff}}}_{\mathit{i}}-\left[\frac{1}{{{{\mathit{r}}_{\mathrm{ff}}}_{\mathit{i}}}^2+{{{\mathit{r}}_{\mathrm{ff}}}_{\mathit{j}}}^2}\frac{\partial }{\partial \lambda }\mathit{U}\right]{{\mathit{r}}_{\mathrm{ff}}}_{\mathit{j}}\right\}\mathit{i}\\ +\left\{\left[\frac{1}{\mathit{r}}\frac{\partial }{\partial \mathit{r}}\mathit{U}+\frac{{{\mathit{r}}_{\mathrm{ff}}}_{\mathit{k}}}{{\mathit{r}}^2\sqrt{{{{\mathit{r}}_{\mathrm{ff}}}_{\mathit{i}}}^2+{{{\mathit{r}}_{\mathrm{ff}}}_{\mathit{j}}}^2}}\frac{\partial }{\partial \varphi }\mathit{U}\right]{{\mathit{r}}_{\mathrm{ff}}}_{\mathit{j}}+\left[\frac{1}{{{{\mathit{r}}_{\mathrm{ff}}}_{\mathit{i}}}^2+{{{\mathit{r}}_{\mathrm{ff}}}_{\mathit{j}}}^2}\frac{\partial }{\partial \lambda }\mathit{U}\right]{{\mathit{r}}_{\mathrm{ff}}}_{\mathit{i}}\right\}\mathit{j}\\ +\left\{\frac{1}{\mathit{r}}\left(\frac{\partial }{\partial \mathit{r}}\mathit{U}\right){{\mathit{r}}_{\mathrm{ff}}}_{\mathit{k}}+\frac{\sqrt{{{{\mathit{r}}_{\mathrm{ff}}}_{\mathit{i}}}^2+{{{\mathit{r}}_{\mathrm{ff}}}_{\mathit{j}}}^2}}{{\mathit{r}}^2}\frac{\partial }{\partial \varphi }\mathit{U}\right\}\mathit{k}\end{array}$$

球座標における次の偏微分が与えられている。

$$\begin{array}{l}\frac{\partial }{\partial \mathit{r}}\mathit{U}=-\frac{\mu }{{\mathit{r}}^2}\sum _{\mathit{l}=2}^{{\mathit{l}}_{\max }}\sum _{\mathit{m}=0}^{\mathit{l}}{\left(\frac{{\mathit{R}}_{\mathrm{cb}}}{\mathit{r}}\right)}^{\mathit{l}}\left(\mathit{l}+1\right){\mathit{P}}_{\mathit{l},\mathit{m}}\left[\sin \left(\varphi \right)\right]\left\{{\mathit{C}}_{\mathit{l},\mathit{m}}\cos \left(\mathit{m}\lambda \right)+{\mathit{S}}_{\mathit{l},\mathit{m}}\sin \left(\mathit{m}\lambda \right)\right\}\\ \frac{\partial }{\partial \varphi }\mathit{U}=\frac{\mu }{\mathit{r}}\sum _{\mathit{l}=2}^{{\mathit{l}}_{\max }}\sum _{\mathit{m}=0}^{\mathit{l}}{\left(\frac{{\mathit{R}}_{\mathrm{cb}}}{\mathit{r}}\right)}^{\mathit{l}}\left\{{\mathit{P}}_{\mathit{l},\mathit{m}+1}\left[\sin \left(\varphi \right)\right]-\left(\mathit{m}\right)\tan \left(\varphi \right){\mathit{P}}_{\mathit{l},\mathit{m}}\left[\sin \left(\varphi \right)\right]\right\}\left\{{\mathit{C}}_{\mathit{l},\mathit{m}}\cos \left(\mathit{m}\lambda \right)+{\mathit{S}}_{\mathit{l},\mathit{m}}\sin \left(\mathit{m}\lambda \right)\right\}\\ \frac{\partial }{\partial \lambda }\mathit{U}=\frac{\mu }{\mathit{r}}\sum _{\mathit{l}=2}^{{\mathit{l}}_{\max }}\sum _{\mathit{m}=0}^{\mathit{l}}{\left(\frac{{\mathit{R}}_{\mathrm{cb}}}{\mathit{r}}\right)}^{\mathit{l}}\left(\mathit{m}\right){\mathit{P}}_{\mathit{l},\mathit{m}}\left[\sin \left(\varphi \right)\right]\left\{{\mathit{S}}_{\mathit{l},\mathit{m}}\cos \left(\mathit{m}\lambda \right)-{\mathit{C}}_{\mathit{l},\mathit{m}}\sin \left(\mathit{m}\lambda \right)\right\}\end{array}$$

${\mathit{P}}_{\mathit{l},\mathit{m}}$ はルジャンドル関数です。

${\mathit{C}}_{\mathit{l},\mathit{m}}$ と ${\mathit{S}}_{\mathit{l},\mathit{m}}$ は正規化されていない調和係数です。

参考文献

[1] ヴァラド、デイビッド。天体力学の基礎と応用、第 4 版。ホーソーン、カリフォルニア州:マイクロコズムプレス、2013年。

[2] Vepa, R.、「自律型宇宙船とロボットのダイナミクスと制御」、University Printing House、ニューヨーク、ニューヨーク、米国、2019年。

[3] Stevens, B. L.、F. L. Lewis、「航空機制御とシミュレーション」第2版、John Wiley & Sons、ホーボーケンNJ、2003年。

[4] ゴットリーブ、R.G.、「高速重力、重力部分、正規化重力、重力勾配トルクおよび磁場：「導出、コードおよびデータ」、NASA 契約者レポート 188243 技術レポート、NASA リンドン B. ジョンソン宇宙センター、テキサス州ヒューストン、1993 年 2 月。

[5] Konopliv, A. S.、S. W. Asmar、E. Carranza、W. L. Sjogen、D. N. Yuan.、「ルナプロスペクターミッションの結果としての最近の重力モデル、イカロス」、第150巻、第1号、pp 1-18、2001年。

[6] Lemoine, F. G.、D. E. Smith、D.D. Rowlands、M.T. Zuber、G. A. Neumann、およびD. S. Chinn、「火星グローバルサーベイヤーによる火星の重力場（GMM-2B）の改良解」、Journal Of Geophysical Research、Vol. 106、No. E10、pp 23359-23376、2001年10月25日。

[7] Seidelmann, P.K., Archinal, B.A., A’hearn, M.F. et al. 地図座標と回転要素に関するIAU/IAGワーキンググループ報告書:2006年。Celestial Mech Dyn Astr 98、155–180 (2007)。

[8] スタンディッシュ、E.M.、「JPL惑星および月暦、DE405 / LE405」、JPL IOM 312.F-98-048、パサデナ、カリフォルニア州、1998年。

参考

Spacecraft Dynamics