(57)【要約】【課題】公知の数値解方式に基づいた種々の計算方法と比べて、計算結果が得られるまでに要する計算時間が格段に短く、また格段に精度のよい計算結果が得られる、宇宙航行機の自己汚染進行過程の計算方法を提案する。【解決手段】この方法は、前記宇宙航行機の一般的定義から成る第１入力パラメータ集合を受取るステップと、前記宇宙航行機の軌道データ、物理量、数値に関する制御パラメータから成る第２入力パラメータ集合を受取るステップと、受取った前記第１入力パラメータ集合及び前記第２入力パラメータ集合のデータに基づいて前記宇宙航行機の自己汚染進行過程を計算する自己汚染進行過程計算ステップであって、前記宇宙航行機の表面から放出された気体分子が堆積する堆積量を算出するための基礎方程式の解析解の算出を前記データ処理装置による一度の計算ステップで行うことで前記宇宙航行機の自己汚染進行過程を計算するようにした、自己汚染進行過程計算ステップと、算出した堆積量を出力するステップとを含むことを特徴とする。【選択図】図１
2015113121.tif 000002
【特許請求の範囲】
【請求項１】
データ処理装置により宇宙航行機の自己汚染進行過程を計算する方法であって、
前記宇宙航行機の一般的定義を含む第１入力パラメータ集合を受取るステップと、
前記宇宙航行機の軌道データ、物理量、数値に関する制御パラメータを含む第２入力パラメータ集合を受取るステップと、
受取った前記第１入力パラメータ集合及び前記第２入力パラメータ集合のデータに基づいて前記宇宙航行機の自己汚染進行過程を計算する自己汚染進行過程計算ステップであって、前記宇宙航行機の表面から放出された気体分子が堆積する堆積量を算出するための基礎方程式の解析解の算出を前記データ処理装置による一度の計算ステップで行うことで前記宇宙航行機の自己汚染進行過程を計算するようにした、自己汚染進行過程計算ステップと、
算出した堆積量を出力するステップと、
を含むことを特徴とする方法。
【請求項２】
前記自己汚染進行過程計算ステップは、前記宇宙航行機を複数の表面の集合体として表し、前記複数の表面の各々を夫々が固有の気体放出特性を有する複数のメッシュエレメントに分割したときの、前記宇宙航行機の所定の１つの表面の第i番メッシュエレメントに対する前記宇宙航行機の所定の１つの表面の第j番メッシュエレメントのビューファクタFjiを算出することを特徴とする請求項１記載の方法。
【請求項３】
前記基礎方程式の解析解の算出は、前記宇宙航行機の前記複数の表面の放出係数及び再放出係数が区分的に定数である場合の堆積量を、下記の基礎方程式を用いて算出し、
【数１】
2015113121.tif 000055
ここで、yは時刻tにおける堆積量であり、Aは前記ビューファクタFjiと前記再放出係数とから算出された定数行列であり、bは前記放出係数とビューファクタと初期堆積量とから構成されたベクトルである、ことを特徴とする請求項２記載の方法。
【請求項４】
前記堆積量の算出は、下記の式に従って積分式の計算を実行することにより、第ｎ番期間の時刻tnにおける堆積量を算出し、
【数２】
2015113121.tif 000056
ここで、yn-1は先行するn-1回の計算ステップにおいて夫々に算出された堆積量の累積値であり、Anは前記ビューファクタFjiと前記再放出係数とから算出された区分的に定数の定数行列であり、b(s)は前記放出係数と前記ビューファクタと時刻sにおける初期堆積量とから構成されたベクトルである、ことを特徴とする請求項３記載の方法。
【請求項５】
前記基礎方程式の解析解の算出は、前記宇宙航行機の前記複数の表面の放出係数及び再放出係数が時間依存性を有する場合の堆積量を、その時間依存性が一次関数であるものとした下記の基礎方程式を用いて算出し、
【数３】
2015113121.tif 000057
ここで、yは時刻tにおける堆積量、A及びBは時間依存性を有する前記ビューファクタFjiと前記再放出係数とから算出された定数行列であり、bは時間依存性を有する前記放出係数と前記ビューファクタと初期堆積量とから構成されたベクトルである、ことを特徴とする請求項２記載の方法。
【請求項６】
前記堆積量の算出は、下記の積分式の計算を実行することにより、前記行列A及びBが定数行列である0≦t≦tAの期間内に含まれる時刻tにおける堆積量を算出し、
【数４】
2015113121.tif 000058
ここで、積分定数cは時刻tにおける堆積量がc=y0であるとした境界条件y(0)=y0を挿入することにより決定された定数であり、b(s)は時間依存性を有する前記放出係数と前記ビューファクタと時刻sにおける初期堆積量とから構成されたベクトルである、ことを特徴とする請求項５記載の方法。
【請求項７】
請求項１乃至６の何れか１項記載の方法を実行するためのコンピュータ・プログラム。
【請求項８】
請求項７記載のコンピュータ・プログラムを格納した記録媒体。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は宇宙航行機の自己汚染進行過程の計算方法に関する。
【背景技術】
【０００２】
人工衛星などの宇宙航行機は、真空下にあってその表面から気体分子が脱離して放出され、特に表面に形成されている被覆層は多くの気体分子を放出し、放出された気体分子による宇宙航行機の自己汚染が進行する。この自己汚染によって表面特性に不都合な変化が生じることがあり、例えば光学表面であれば、そのスペクトル特性や方向性特性が劣化するおそれがあり、また非光学表面であれば、その放出特性、吸収特性、反射特性などが劣化するおそれがある。非特許文献１には、ロゼッタ彗星探査機ミッションの遂行中に測定された気体分子放出及び汚染と、それによって発生した様々な問題とが論じられている。
【０００３】
そのような不都合な表面特性の変化は、場合によってはスペース・ミッションにとって致命的なものとなりかねず、そのため重要なスペース・ミッションでは、そのミッションの開始前に、自己汚染進行過程を計算することで不都合な表面特性の変化を把握するようにしており、その計算を行うための、数値計算法に基づいた様々な解析モデルが開発されている。現代の宇宙航行機には様々な材料物質が用いられており、自己汚染の進行速度は汚染を引き起こす物質ごとに大きく異なっている。しかるに、公知の解析モデルに採用されている自己汚染進行過程を計算するための種々の数値計算法はいずれも長い計算時間を要するものであり、特に自己汚染の進行速度が速い場合にはその計算時間が非常に長くなる。更に加えて、それら数値計算法では多くの場合、かなりの計算誤差を生じる。
【０００４】
欧州宇宙機関（ＥＳＡ）は、放出された気体分子による宇宙航行機の表面汚染のモデルを生成するためのソフトウェアを提供しており、このソフトウェアは、その名称を「ＣＯＭＯＶＡ」（COntamination Modelling Outgassing & Vent Analyses tool）という（同ソフトウェアの詳細はインターネット上のサイトhttp://www-mip.onera.fr/comova/に開示されている）。ＣＯＭＯＶＡを使用することで宇宙航行機の表面汚染を計算することができるが、ただしその計算には長時間を要する。ＣＯＭＯＶＡは更に、「ＥＳＡＢＡＳＥ２宇宙工学解析ツール」という名称のソフトウェア（同ソフトウェアの詳細はインターネット上のサイトhttp://esabase2.net/product/esabase2-comova/に開示されている）に組込んで使用することで表面に堆積する汚染物質の堆積量を計算することができ、その堆積量の計算には、ビューファクタを求めるための計算と、表面温度の影響の計算とが含まれる。
【先行技術文献】
【非特許文献】
【０００５】
【非特許文献１】Ｍ．ハシック他著、「宇宙航行機における気体分子放出−過小評価されている現象」（“Spacecraft outgassing, a largely underestimated phenomen”, M. Hassig et al., 2nd International Conference on Space Technology (ICST), 15-17 Sept. 2011, pages 1-4 Print ISBN: 978-1-4577-1874-8）
【発明の概要】
【発明が解決しようとする課題】
【０００６】
本発明の目的は、公知の数値解方式に基づいた種々の計算方法と比べて、計算結果が得られるまでに要する計算時間が格段に短く、また格段に精度のよい計算結果が得られる、宇宙航行機の自己汚染進行過程の計算方法を提案することにある。
【課題を解決するための手段】
【０００７】
上記目的は独立請求項の主題により達成される。尚、従属請求項は実施の形態を示したものである。
【０００８】
本発明は、宇宙航行機の自己汚染進行過程を計算するための新規なアルゴリズムを提案するものであり、このアルゴリズムは、宇宙航行機の様々な材料物質の気体分子の放出及び自己汚染という問題の根本的性質を表す数式の解析解に基づく。
【０００９】
本発明の１つの実施の形態は、データ処理装置により宇宙航行機の自己汚染進行過程を計算する方法に関する。この方法は、前記宇宙航行機の一般的定義から成る第１入力パラメータ集合を受取るステップと、前記宇宙航行機の軌道データ、物理量、数値に関する制御パラメータから成る第２入力パラメータ集合を受取るステップと、受取った前記第１入力パラメータ集合及び前記第２入力パラメータ集合のデータに基づいて前記宇宙航行機の自己汚染進行過程を計算する自己汚染進行過程計算ステップであって、前記宇宙航行機の表面から放出された気体分子が堆積する堆積量を算出するための基礎方程式の解析解の算出を前記データ処理装置による一度の計算ステップで行うことで前記宇宙航行機の自己汚染進行過程を計算するようにした、自己汚染進行過程計算ステップと、算出した堆積量を出力するステップとを含むことを特徴とする。表面特性に応じて、即ち、例えば、放出係数及び再放出係数が区分的に定数であるかそれとも時間依存性を有しているかに応じて、前記宇宙航行機の表面から放出された気体分子が堆積する堆積量を算出するため前記基礎方程式として、適当な基礎方程式を選択することができ、また、その選択した基礎方程式の解析解の算出を一度の計算ステップで行うことができるため、前記基礎方程式の解析解を算出する上で、時間のかかる無用な中間ステップを反復実行する必要がない。
【００１０】
前記自己汚染進行過程計算ステップにおいては、前記宇宙航行機を複数の表面の集合体として表し、それら表面の各々を夫々が固有の気体放出特性を有する複数のメッシュエレメント（以下「ノード」と呼ぶ）に分割したときの、前記宇宙航行機のある１つの表面の第i番メッシュエレメントに対する前記宇宙航行機のある１つの表面の第j番メッシュエレメントのビューファクタFjiを算出するようにするとよい。
【００１１】
前記基礎方程式の解析解の算出においては、前記宇宙航行機の前記複数の表面の放出係数及び再放出係数が区分的に定数である場合の堆積量を、下記の基礎方程式を用いて算出するようにするのもよく、
【数１】
2015113121.tif 000003
ここで、yは時刻tにおける堆積量であり、Aは前記ビューファクタFjiと前記再放出係数とから算出された定数行列であり、bは前記放出係数と前記ビューファクタと初期堆積量とから構成されたベクトルである。またその場合に、前記堆積量の算出においては、下記の式に従って積分式の計算を実行することにより、第ｎ番堆積期間の時刻tnにおける堆積量を算出するようにするのもよく、
【数２】
2015113121.tif 000004
ここで、yn-1は先行するn-1回の計算ステップにおいて夫々に算出された堆積量の累積値であり、Anは前記ビューファクタFjiと前記再放出係数とから算出された区分的に定数の定数行列であり、b(s)は前記放出係数と前記ビューファクタと時刻sにおける初期堆積量とから構成されたベクトルである。
【００１２】
或いはまた、前記基礎方程式の解析解の算出においては、前記宇宙航行機の前記複数の表面の放出係数及び再放出係数が時間依存性を有する場合の堆積量を、その時間依存性が一次関数であるものとした下記の基礎方程式を用いて算出するようにするのもよく、
【数３】
2015113121.tif 000005
ここで、yは時刻tにおける堆積量、A及びBは時間依存性を有する前記ビューファクタFjiと前記再放出係数とから算出された定数行列であり、bは時間依存性を有する前記放出係数と前記ビューファクタと初期堆積量とから構成されたベクトルである。またその場合に、前記堆積量の算出においては、下記の積分式の計算を実行することにより、前記行列A及びBが定数行列である0≦t≦tAの期間内に含まれる時刻tにおける堆積量を算出するようにするのもよく、
【数４】
2015113121.tif 000006
ここで、積分定数cは時刻t=0における堆積量がc=y0であるとした境界条件y(0)=y0を挿入することにより決定された定数であり、b(s)は時間依存性を有する前記放出係数と前記ビューファクタと時刻sにおける初期堆積量とから構成されたベクトルである。
【００１３】
本発明の更なる実施の形態はコンピュータ・プログラムに関するものであり、このコンピュータ・プログラムは、本明細書に記載する本発明に係る宇宙航行機の自己汚染進行過程の計算方法を実行するものであって、コンピュータにより実行されることで本発明に係る宇宙航行機の自己汚染進行過程を計算することを可能にするものである。かかるコンピュータ・プログラムは、宇宙航行機に関する計算及びシミュレーションを実行するためのソフトウェア・スイートに組込むのに適したものである。
【００１４】
本発明の更なる実施の形態によれば、本発明に係るコンピュータ・プログラムを格納した記録媒体が提供される。この記録媒体は、例えばＣＤ−ＲＯＭ、ＤＶＤ、メモリーカード、ディスク媒体、及びその他のデータ格納媒体であって、電子的にアクセス可能にコンピュータ・プログラムを格納するのに適した記録媒体である。
【００１５】
本発明の以上の態様は、以下に記載する実施の形態についての説明によって更に明らかとなり、またその他の態様も明らかとなる。
【００１６】
以下に本発明を幾つかの具体的な実施の形態に即して更に詳細に説明する。ただし本発明はそれら具体的な実施の形態のみに限定されるものではない。
【図面の簡単な説明】
【００１７】
【図１】２つの表面i及びjによる自己汚染の進行過程を示した模式図である。
【図２】２つのグラフを示した図であり、一方のグラフは第１ノード及び第２ノードにおける堆積量を示し、他方のグラフは合計堆積量を正規化して示しており、それら堆積量は、放出係数及び再放出係数が共に定数の２枚の板材から成る試験例に関して本発明に係る方法を用いて算出した堆積量である。
【図３】合計堆積量の算出値における誤差を特異値の個数の関数として示したグラフである。
【図４】連続する２回の計算ステップで夫々に算出された合計堆積量の算出値の間の誤差を特異値の個数の関数として示したグラフである。
【発明を実施するための形態】
【００１８】
以下の説明では、機能的に同等ないし同一の要素を示すために同一の参照符号を用いることがある。また、以下の説明において具体的な数値を提示することがあるが、その数値は具体例を示すものに過ぎず、本発明を限定するものではない。
【００１９】
（０）宇宙航行機の自己汚染進行過程の計算の基礎
宇宙航行機の自己汚染進行過程を計算するために、また特に、複数種類の材料物質が使用されている宇宙航行機の自己汚染進行過程を計算するために、その宇宙航行機を複数の表面の集合体として表し、更に、それら表面の各々を、夫々が固有の気体放出特性を有する複数のノード（メッシュエレメント）に分割する。それら複数のノードは、夫々が汚染物質を放出し、またノードどうしの間で汚染物質のやり取りを行う。宇宙航行機のある表面から（即ち使用材料から）放出された気体分子は、宇宙航行機の全ての表面へ向かって移動し、到達先の表面を汚染する汚染物質となる。
【００２０】
データ処理装置により宇宙航行機の自己汚染進行過程を計算するというとき、そのデータ処理装置とは、例えば、本発明に係る方法の基本ステップを実施するコンピュータ・プログラムを実行するコンピュータなどである。データ処理装置により宇宙航行機の自己汚染進行過程を計算するために、データ処理装置は２組の入力パラメータ集合を受取る。それらのうち、第１入力パラメータ集合は、宇宙航行機の一般的定義を含むものであり、この一般的定義には、例えば、宇宙航行機の形状や宇宙航行機の表面の特性などが含まれ、また特に、宇宙航行機の表面の（表面に被覆層が形成されている場合にはその被覆層の）気体分子の放出特性及び再放出特性や温度などが含まれる。また、第２入力パラメータ集合は、宇宙航行機の軌道データ、物理量、数値に関する制御パラメータを含むものである。宇宙航行機の軌道データの典型的なものとしては、軌道定義データ、軌道定義データに内挿を施すための内挿データ、それに運動学的データなどがある。物理量に関する制御パラメータとしては、大気定義データや輸送定義データなどがある。数値に関する制御パラメータとしては、レイトレーシング及びビューファクタに関するパラメータ、汚染及び気体放出により影響を受ける表面に関するパラメータ、及び／または、ボリュームメッシュデータなどがある。
【００２１】
データ処理装置は、宇宙航行機の自己汚染進行過程を計算するための入力データを受取ったならば、その受取った入力データに計算処理を施す。その計算処理には、宇宙航行機の表面から放出された気体分子が堆積する堆積量を算出するための基礎方程式の解析解の算出が含まれる。また、その計算処理は次の２つの基本ステップに分けられる。
１．ビューファクタを算出するステップ。ビューファクタとは、第i番ノードから放出される物質の全質量のうち第j番ノードに到達する物質の質量の割合である。
２．堆積量を表す積分方程式の解を算出するステップ。この積分方程式の解の算出は基礎方程式の解析解を算出するという形で行う。
【００２２】
第１のステップでは、第i番表面メッシュエレメント（ノード）と第j番表面メッシュエレメント（ノード）との間のビューファクタFijを算出する。ビューファクタFijをマクロのレベルで定義するならば、第i番表面メッシュエレメント（ノード）から放出されて第j番表面メッシュエレメントに到達した物質の質量（または分子数）を、第i番表面メッシュエレメント（ノード）から放出された物質の全質量（または全分子数）で除した値であると定義することができる。一方、ミクロのレベルで定義するならば、ビューファクタFijは、第i番ノードから放出された分子が第j番ノードに到達する確率であると定義することができる。
【００２３】
図１に示すように、第i番ノードから第j番ノードへ移転する物質の質量は様々な質量束で表すことができ、それら質量束のうちには、直接的に移転する質量の質量束（直接質量束）、他の汚染物質分子により散乱されて移転する質量の質量束（自己散乱質量束）、環境中の気体の流束により散乱されて移転する質量の質量束（環境気体散乱質量束）、それに、宇宙航行機の表面から再放出されて移転する質量の質量束（再放出質量束）が含まれる。
【００２４】
ビューファクタFijの値を算出するにはモンテカルロ・レイトレーシング法を用いるとよい。それについて手短に説明すると、先ず、放出される個々の気体分子をレイで表す。そうすると、第i番ノードは多数のレイをランダムに射出していることになる。このレイの射出に関しては、射出点のサンプリングを等方分布則に従って行い、射出方向のサンプリングをランバート則に従って行う。また、平均自由行程λaを、剛体球モデルを用いて断面積を算出するロバートソンの式に従って算出する。更に、各々のレイについての衝突距離dcのサンプリングを衝突確率に従ってランダムに行い、その衝突確率P(x)は下記の式で表される。
【数５】
2015113121.tif 000007
この式において、xは距離であり、λaは平均自由行程である。衝突距離dcとノード間距離djとの比較によって、当該分子が散乱されるかそれとも第j番ノードに到達するかが決定される。dc＜djであるならば当該分子は衝突し、そうでなければ当該分子は衝突しない。衝突した分子はその衝突点からランダムな方向に再放出され、そして第k番ノードに到達してFikに寄与することになる。
【００２５】
ビューファクタFjiを算出したならば、次に、堆積量（ノードに堆積した物質の質量）を算出する。それについて手短に説明すると、以下の通りである。第i番ノード、第j番ノード、第k番ノード、…、によってネットワークが形成されており、そのネットワークに含まれるノードのうちの２つのノードから成る各ペア(i,j)、(i,k)、(j,k)、…、は、そのペアの２つのノードの関連度が、対応するビューファクタFij、Fik、Fjk、…、により表されている。また、各ノードは、汚染物質を放出し、汚染物質が堆積し、堆積した汚染物質を再放出する。従って、堆積量を求める問題は一階微分方程式系となり、本発明の解法によってこの一階微分方程式系の解を算出することが可能であり、それについて以下に詳細に説明する。
【００２６】
以下に説明する本発明の解法における計算ステップは、例えば上で言及したＣＯＭＯＶＡソフトウェアの追加モジュールないしプラグインなどの、コンピュータ・プログラムに組込むのにも適しており、また、夫々に気体の放出特性ないし脱離特性が異なる種々の材料物質が使用されている宇宙航行機の自己汚染進行過程を計算するために必要とされる計算時間を、大幅に短縮し得るものである。尚、かかる計算ステップの背景と幾つかの方程式の解法についての説明を、実施形態の説明の末尾に「付記」として記すことにする。
【００２７】
（１）放出係数及び再放出係数が区分的に定数である場合の堆積量の算出
（１．１）堆積量を算出するための基礎方程式
第i番ノードの堆積量を算出するためにその解を求める方程式は、下記の式（１）で表される（この方程式は「ＯＵＴＧＡＳＳＩＮＧソフトウェア・ユーザーズ・マニュアル」（"OUTGASSING Software User's Manual", internal document identification MOS.UM.93694.ASTR）の第１３６頁に式（７２）として掲載されているものであり、同文献の記載内容はこの言及をもって本願開示に組込まれたものとする。また、この方程式を簡略化した方程式がVan Eesbeek, M.、Zwaal A.共著、「滞留時間に基づく気体放出及び汚染モデル」（"Outgassing and Contamination Model Based on Residence Time", Proceedings of the Third European Symposium on Spacecraft Materials in Space Environment, Noordwijk, Oct. 1985, (ESA SP-232））に記載されている。
【数６】
2015113121.tif 000008
この式（１）において、Fjiは第j番ノード（第１表面メッシュエレメント）から第i番ノード（第２表面メッシュエレメント）へのビューファクタであり、τr,jは第j番ノードの再放出係数であり、τe,jは第j番ノードの放出係数である。また、この式（１）において、係数行列Aijは定数行列である。
【００２８】
（１．２）係数行列が区分的に定数行列である場合の解析解の一般形
理解を容易にするために、式（１）を簡略化した方程式である、下記の式（２）の解の求め方について説明する。
【数７】
2015113121.tif 000009
この式（２）からは、下記の式（３）が得られ、
【数８】
2015113121.tif 000010
なぜならば、
【数９】
2015113121.tif 000011
だからである。
またここで、時刻tがtA≦t≦tBである期間内において、行列Aは定数行列であるものとする。
【００２９】
すると、式（３）から、下記の式が得られる。
【数１０】
2015113121.tif 000012
この式を、時刻tAから時刻tまで積分することにより、下記の式（４）が得られる。
【数１１】
2015113121.tif 000013
この式（４）に、時刻t＝tAにおける堆積量がy(tA) = yAであるという境界条件を挿入することにより、積分定数cの値を求めることができ、その値は下記の式（５）で表される。
【数１２】
2015113121.tif 000014
従って、完全な解は下記の式（６）の形を取る。
【数１３】
2015113121.tif 000015
【００３０】
ここで、係数行列Aが、連続する個々の期間において夫々に区分的に定数行列であれば、以上の解析プロセスを連続する個々の期間において夫々に実行することができる。そこで、
0≦t1≦t2、且つ
A(t) ≡ A1＝定数行列  ただし0≦t≦t1、且つ
A(t) ≡ A2＝定数行列  ただしt1＜t≦t2、とする。
【００３１】
こうしてt1＜t≦t2 の期間内の堆積量を算出すると、その堆積量は下記の式（７）で表される。
【数１４】
2015113121.tif 000016
この式（７）の中のy1は、下記の式で表される。
【数１５】
2015113121.tif 000017
この式の中のy0は、時刻t＝0における堆積量であり、通常は0である。
【００３２】
以上の式を、連続するn個の期間にまで一般化することにより、第n番期間に含まれる時刻tnにおける堆積量が得られ、その堆積量は下記の式（８）で表される。
【数１６】
2015113121.tif 000018
この式（８）において、yn-1は、先行するn-1回の計算ステップの各々において算出した堆積量を合計した累積堆積量である。
【００３３】
（１．３）係数行列が区分的に定数行列である場合の解析解の明示形
式（８）の積分を実行するためには、式（２）の固有値分解を用いる必要があり、その固有値分解は下記の式（９）で表される。
【数１７】
2015113121.tif 000019
この式（９）において、行列Pは定数行列Aの（複数の）固有ベクトルを成分とする行列であり、λmは定数行列Aの（複数の）固有値である。ここで、b(s) を式（１）の生成項とするならば、下記の式（１０）に示すように積分が可能である。
【数１８】
2015113121.tif 000020
【００３４】
この式（１０）は、時刻が0からtまでの第１番期間における期間内堆積量の時間依存性を表している。この式（１０）において、W0d,mは第m番ノードの時刻t＝0における堆積量であり、通常は0である。tn-1≦t≦tnで表されるその他の夫々の期間についても、それら期間の期間内堆積量を表す式は、この式（１０）と同じ形を取り、従って、下記の式（１１）で表される。
【数１９】
2015113121.tif 000021
【００３５】
この式（１１）では、期間tn-1≦t≦tnの期間内堆積量が、期間0≦t≦(tn-tn-1）の期間内堆積量にマッピングされており、この式（１１）において、Wnd,mは時刻tn-1までの累積堆積量である。また、Wn-1e,lは修正放出強度であり、下記の式（１２）で表される。
【数２０】
2015113121.tif 000022
また、λ(n)jと、τ(n)e,lとは、期間tn-1≦t≦tnの期間内温度から算出される。
【００３６】
（１．４）特異点解析
式（１０）は特異点を含むものと考えられる。なぜならば、パラメータの値の組合せによっては、λjτe,l＋1の値が0になることがあり得るからである。そこで、そのような特異点が除去可能であることを以下に証明する。
【００３７】
下記の式（１３）が成り立つことは明らかである。
【数２１】
2015113121.tif 000023
【００３８】
ここで、任意のtの値（ただし、その値は108秒＝６年を超えないものとする）に対してλj＋(1/τe,l)≪tが成り立つならば、式（１３）は展開可能であり、展開することで下記の式（１４）が得られ、
【数２２】
2015113121.tif 000024
この式（１４）の右辺は有限である。
【００３９】
（１．５）試験例
（１．５．１）放出係数及び再放出係数が定数である２枚の板材に関する試験
相互に完全なビジビリティを有する２枚の板材の間で物質が移転する場合について式（１）の解を求める試験を行った。この試験では、以上に説明した解析解を利用して、マスワークス社から提供されているＭＡＴＬＡＢ（登録商標）というコンピュータ・プログラムが、ある期間に亘って変化する式（１）で示された堆積量を算出し得るか否かを判定した。その期間を表す時間軸は乗用対数目盛で２０桁以上の長さに及ぶものとした。また、その常用対数目盛の時間軸上に等間隔で位置する時刻の夫々における堆積量を算出させ、それによって、宇宙航行機の供用開始直後に堆積量のピークが発生した場合であってもユーザがそのピークを一目で視認でき、且つ、宇宙航行機の寿命終了時点での堆積量の値もユーザが容易に読み取れるようにした。
【００４０】
使用したパラメータは以下の通りである。
【数２３】
2015113121.tif 000025
尚、上記パラメータのうち、N＝100としているのは、計算ステップを100回実行したこと（時間軸の全長を100等分する時刻の各々で計算ステップを実行したこと）を表している。
【００４１】
図２のグラフに示したように、算出された堆積量は滑らかに変化しており、時間軸の末端付近でも十分に精細な解像度が得られている。また、２枚の板材の相対位置は、物質が逸失しないような相対位置に設定したため、時間軸の末端に至っても物質はその全質量が保存されている。
【００４２】
堆積量の算出値は、解析解を求める式を最大桁の１６桁にした場合も同一であった。
【００４３】
（２）放出係数及び再放出係数が区分的に一次関数的に変化する場合の堆積量の計算
（２．１）放出係数が時間依存性を有する場合の放出則
放出係数が時間依存性を有さない場合（定数である場合）には、使用材料から放出される物質の放出量は下記の式（１５）で表される。
【数２４】
2015113121.tif 000026
この式（１５）において、W0eは放出可能な物質の全質量である。
【００４４】
また、この式（１５）は、一次反応速度式である下記の式（１６）の解である。
【数２５】
2015113121.tif 000027
この式（１６）において、τeは時間依存性を有さない放出係数であり、アレニウス則に従って下記の式（１７）により決定されるものである。
【数２６】
2015113121.tif 000028
【００４５】
以下の説明では、時刻tがtA≦t≦tBである期間内において、温度Tが一次関数的に変化するものとする。これに関しては、欧州宇宙技術研究センター（ＥＳＴＥＣ）及びフランス国立宇宙研究センター（ＣＮＥＳ）が提唱している放出則でも、またアレニウス則でも、温度Tは指数関数的に変化するものとしている。しかしながら、tA≦t≦tBの期間内に限定するならば、式（１７）を下記の式（１８）で近似することができ、それら２つの式の差分は有限である。
【数２７】
2015113121.tif 000029
【００４６】
この近似式（１８）は、ＥＳＴＥＣ及びＣＮＥＳが提唱しているモデルにも適用可能であり、適用する際には、有限の差分商を算出するときと同様に、用いようとする放出則モデルに単に時刻tA及びtBにおける温度の値を設定するだけでよい。
【００４７】
式（１８）を式（１６）に代入して積分を行うことにより、下記の式（１９）が得られる。
【数２８】
2015113121.tif 000030
この式（１９）の右辺を展開することによって、τe・（式中にτeの上部にドットを付して表したτe）の１階導関数即ち傾きで表される下記の式（２０）が得られる。
【数２９】
2015113121.tif 000031
この式（２０）の右辺は、τeの傾きが小さくなるほど旧来の放出則に近付いて行き、その傾きがゼロなったならば旧来の放出則に還元される。
【００４８】
（２．２）堆積量を算出するための基礎方程式
放出係数のみならず、再放出係数も同様に、時間依存性を有するものと考える必要がある。これは、式（１）の中の行列Aijを、時間依存性を有する行列として取扱うことを意味する。その時間依存性は、式（１８）を用いることにより、tA≦t≦tBの期間内における時刻tの一次関数で表すことができ、それによって下記の式（２１）が得られる。
【数３０】
2015113121.tif 000032
【００４９】
この式（２１）も式（１８）と同様に、用いようとする放出則モデル（例えばＣＮＥＳないしＥＳＴＥＣが提唱しているモデルなど）に対する依存性を有しておらず、なぜならば、この式（２１）の中の再放出係数が関わるのは、時刻tA及びtBにおける温度の値だけだからである。この式（２１）に前章で導出した式を代入することによって、堆積量を算出するための方程式が下記の式（２２）の形で得られる。
【数３１】
2015113121.tif 000033
【００５０】
（２．３）近似解の一般形
上記の式（２２）は、時間依存性を有する係数を含む非均質常微分方程式系である。先の「放出係数及び再放出係数が区分的に定数である場合の堆積量の算出」の章で説明したのと同様にして、係数を変化させる方法を適用するためには、均質系の基本解系と呼ばれる、全ての解を含む解ベクトルが必要である。しかしながら通常は、そのような基本解系を解析的に算出することは不可能である。そこで、基本解系を用いるのではなく、この常微分方程式系を先ず、積分方程式の形に変え、続いて、台形則を用いてその積分方程式の近似解を求めるようにする。
【００５１】
式（２）と同様の表記法を用いることによって、式（２２）を下記の式（２３）のように表すことができる。
【数３２】
2015113121.tif 000034
この式（２３）から、下記の式（２４）が得られる。
【数３３】
2015113121.tif 000035
またここで、時刻tが0≦t≦tAである期間内において、行列A及びBは定数行列であるものとする。
すると、式（２４）から、下記の式（２５）が得られる。
【数３４】
2015113121.tif 000036
この式（２５）に、時刻t＝0における堆積量がy(0)＝y0であるという境界条件を挿入することによって積分定数cの値が得られ、その値は下記の式（２６）で表される。
【数３５】
2015113121.tif 000037
【００５２】
（２．４）近似解の明示形
上記の式（２５）は、下記の式（２７）のように書き換えられる。
【数３６】
2015113121.tif 000038
更に、行列Aを対角化する。それを添字記法で表すと、下記の式（２８）となる。
【数３７】
2015113121.tif 000039
【００５３】
この式（２８）の右辺の第２項に、台形則による近似を適用すると、同項は、下記の式（２９）で表される。
【数３８】
2015113121.tif 000040
この近似の適用により、下記の式（３０）が得られる。
【数３９】
2015113121.tif 000041
式（３０）をyについて解くことで、下記の式（３１）が得られる。
【数４０】
2015113121.tif 000042
【００５４】
ただし、この式（３１）の中の逆行列を計算によって明示形の行列として求めるようなことは、通常は行わない。それよりも、ガウスの消去法を用いて式（３０）の解を求める方が、計算処理が短時間で完了し、精度も良好となる。
【００５５】
そこで、先ず、式（３０）の中のyci(t)を求めるための計算処理について説明する。式（３０）が導出された時点では、もはや放出則は指数関数を含んでいないため、積分方程式の解を解析的に求めることはできない。ただし、その替わりに、下記の式（３２）が得られている。
【数４１】
2015113121.tif 000043
【００５６】
この式（３２）に従って、各々の添字ペア(j,n)ごとに積分値を算出する。そして、算出した積分値に、先ず、下記の項を乗じる。
【数４２】
2015113121.tif 000044
更に続いて、行列P-1jmFTmnの対応する成分を乗じる。尚、積分値を評価するには、不完全超幾何関数の計算を行わねばならず、それには長時間を要することがある。
【００５７】
（２．５）試験例
（２．５．１）放出係数及び再放出係数が時間依存性を有する２枚の板材に関する試験
下記のパラメータを用いて、式（３１）及び式（３２）の解を算出する試験を行った。
【数４３】
2015113121.tif 000045
【００５８】
下記の表に、tの値を様々に設定して得られた試験結果を示した。括弧内の一対の値は左側が式（３１）の解を算出するのに要した時間、右側が式（３２）の解を算出するのに要した時間である。更に、ＭＡＴＬＡＢ（登録商標）の「ＯＤＥ (Ordinary Differential Equation: 常微分方程式)４５ソルバ」及び「ＯＤＥ１５ｓソルバ」を用いて解の検証を行った。
【表１】
2015113121.tif 000046
【００５９】
（３）近似的不完全行列の取扱い
（３．１）擬似逆行列の算出
行列Aは一般的に、非対称行列である上に極めて大型の行列となるため、対角化が困難なことがある。ただし、成分数が数桁以上にも及ぶ大型行列では、その行列の固有ベクトルどうしが互いに非常に近似していることがある。
【００６０】
実際には、ある行列が厳密に数学的な意味において不完全行列となる可能性は極めて低い（これについては末尾の「付記」を参照されたい）のであるが、しかしながら、その行列の固有ベクトルどうしが互いにあまりにも近似しているために、それら固有ベクトルが互いに異なっていると判断するのが困難なことが多々ある。そのような行列Pでは、条件数が膨大となることから逆行列を求めることができず、また、丸め誤差が勝手に拡大することから算出結果は信頼性の乏しいものとなる。ＭＡＴＬＡＢ（登録商標）にはrank関数が装備されており、このrank関数によって、行列に含まれる線形独立な列の個数及び線形独立な行の個数を算出することができる。
【００６１】
ある行列Pのrankとは、その行列Pの特異値のうち、その値が下記の式（３３）で表されるtolの値より大きい全ての特異値の個数であると定義される。
【数４４】
2015113121.tif 000047
この式（３３）において、Nは行列Pの次元数であり、siiは行列Pの特異値を成分とする対角行列である。
【００６２】
行列Pのrankをrank(P)で表すとき、rank(P)＜N であるならば、複素行列Pの逆行列を求めることはできないが、ただし、その逆行列を擬似逆行列P+によって近似することが可能である。
【００６３】
ここで、行列Pの特異値分解が、下記の式（３４）で与えられるならば、
【数４５】
2015113121.tif 000048
擬似逆行列P+は、下記の式（３５）で示したように構成することができる。
【数４６】
2015113121.tif 000049
この式（３５）において、S-1truncatedは、行列Sのトランケーテッド逆行列であり、これは、許容誤差の値よりも小さな全ての特異値を0に設定した上で、残余の特異値を成分とする行列の逆行列を求めたものである。
【００６４】
（３．２）特異値の個数が堆積量の算出値に及ぼす影響
行列Pの真正の逆行列に代えて式（３５）に示した擬似逆行列P+を用いることにより生じる誤差の大きさを調べるために、実行中の解析から取出したrank数がN=2161の非不完全問題行列Aijを用いて誤差についての解析を行った。この行列は対角化が可能であったため、様々な擬似逆行列を用いて得られる堆積量の算出値を検証するための基準値を得ることができた。計算ステップを21回実行して、各回の計算ステップにおいて全てのノードの堆積量を算出した。計算ステップを１回実行する度に、擬似逆行列を構成するために用いる特異値を100個ずつ増やして行った。22回目に実行した計算ステップでは、全ての特異値を用いることになった。この22回目に実行した計算ステップで算出した堆積量の算出値を基準値とし、その基準値を行列Pの真正の逆行列inv(P) を用いて算出した堆積量の算出値と比較したところ、相対誤差は1e-10以下であった。
【００６５】
以下に説明するグラフは、シミュレーション終了時点での全てのノードの堆積量を合計した合計堆積量における相対誤差の大きさを、擬似逆行列を構成するために用いた特異値の個数の関数として表したものである。
【００６６】
図３に示したグラフは、合計堆積量の算出値における誤差を特異値の個数の関数として示したグラフである。
【００６７】
図３のグラフから明らかなように、2161個のうちの2100個の特異値を用いたときの相対誤差は1e-9以下であった。
【００６８】
問題行列が実際に不完全行列である場合もあるため、それに備えて収束性解析を行っておく必要がある。この収束性解析は、より多くの特異値を用いて得られた算出値と、より少ない特異値を用いて得られた算出値とを比較することにより行う。ある回の計算ステップの実行で得られた算出値と、それに続く次の回の計算ステップの実行で得られた算出値との間の相対誤差が、所定のスレショルド値を下回ったならば、そのことによって、計算ステップで用いた擬似逆行列が適切なものであったことが分かる。以下に説明するグラフは、計算ステップを１回実行する度に、次回の計算ステップの実行に用いる擬似逆行列を構成する特異値の個数を100個ずつ増やして行ったときに、連続する２回の計算ステップの実行により夫々に算出された合計堆積量の間の相対誤差が、どのように変化して行くかを示したグラフである。このグラフの右端（x軸＝2100）の値からは、2000個の特異値を用いたときに算出された合計堆積量に対する、2100個の特異値を用いたときに算出された合計堆積量の相対誤差が、2e-5以下であることが分かる。
【００６９】
図４に示したグラフは、連続する２回の計算ステップで夫々に算出された合計堆積量の算出値の間の誤差を特異値の個数の関数として示したグラフである。
【００７０】
（４）付記
（４．１）行列の対角化
ユセフ・サード著、「大規模固有値問題のための数値計算方法」（Yousef Saad, "Numerical methods for large eigenvalue problems", Cambridge University Press, 2nd ed.）に、以下の記載がある。
１．行列Aのある固有値λが特性多項式の重複度μの根であるとき、その固有値λは代数的重複度μを有する。
２．ある固有値が代数的重複度に関わる固有値であるならば、その固有値は単純固有値である。非単純固有値は複数である。
３．行列Aのある固有値λに対応した独立固有ベクトルの最大数がγであるならば、その固有値λは幾何的重複度γを有する。
４．ある行列の少なくとも１つの固有値の代数的重複度が２以上であるならば、その行列は退化行列である。
５．ある固有値の代数的重複度がその固有値の幾何的重複度に等しいならば、その固有値は半単純固有値である。半単純固有値ではない固有値は不完全固有値である。
以下の２つの立言は同等である。
・ｎ次元のある行列が対角化可能であるのは、その行列の線形独立性を有する固有ベクトルの個数がｎ個であるときに限られる。
・ある行列が対角化可能であるのは、その行列の全ての固有値が半単純固有値である場合に限られる。
従って、以下の３つの場合の判別が可能である。
１．全ての固有値が互いに異なる場合：この場合に行列は対角化可能である。
２．全ての固有値が半単純固有値であり、従って、代数的重複度μ＞１であるが、複数の固有値の各々がμ個の線形独立性を有する固有ベクトルを有している場合：この場合に行列は対角化可能である。
３．少なくとも１個の固有値が半単純固有値ではない場合：この場合には行列は対角化可能ではなく、不完全行列となる。
【００７１】
（４．２）相対マシン精度の定義
ＩＥＥＥ規格における相対マシン精度の定義によれば、所与の倍精度数xに対して、それとは絶対値の異なる、それに続く次の倍精度数x+εは、必ず下式に従う。
【数４７】
2015113121.tif 000050
また、ＭＡＴＬＡＢ（登録商標）では、任意の倍精度数Xに対して、そのXと同一精度を有する、それに続く次の大きさの浮動小数点数をε(X)と定義している。
【００７２】
（４．３）逆行列の精度を向上させる幾つかの方法
本発明に係る解法の構造を理解するためには、逆行列P-1lmについての知識が必要である。ＭＡＴＬＡＢ（登録商標）に組込まれている逆行列生成ルーチンを使用したときに、過大と見られる残差が発生することがあり、特に、問題行列Aが対称行列でない場合に大きな残差が発生しがちであり、残差が過大であると、逆行列を求めるための数値計算が速やかに収束しないことがあり、また、十分な精度が得られないことがある。
【００７３】
（左固有ベクトルの利用して精度を向上させる）
ある問題行列Aの左固有ベクトルは、その問題行列Aの転置行列の右固有ベクトルの転置ベクトルと同一であることを利用して精度を向上させることができ、また、複数の行の夫々が左固有ベクトルである行列は複数の列の夫々が右固有ベクトルである行列の逆行列であることを利用して精度を向上させることができる。
【００７４】
（ＭＡＴＬＡＢ（登録商標）で左固有ベクトルを利用するには）
[Vr,Dr]=eig(A) を実行することである問題行列Aの右固有ベクトルを生成できるのであれば、[Vl,Dl]=eig(A') を実行することでその問題行列Aの転置行列の右固有ベクトルを生成することができ、従って、その問題行列Aの左固有ベクトルの転置ベクトルを生成することができる。また、Dr=Dlが成り立つ。ただし、ＭＡＴＬＡＢ（登録商標）が一連の左固有ベクトルと一連の右固有ベクトルとを同一順序で生成するようにしておく必要がある。そして、Vr*Vl'=diagonalの計算を実行することによって、対角性を検証することができる。尚、単位対角行列を得るためには、固有ベクトルに更にスケーリングを施さねばならないことがある。
【００７５】
（反復実行により精度を向上させる）
ＭＡＴＬＡＢ（登録商標）を用いて行列Pの近似逆行列P~-1を算出したときに、近似逆行列P~-1と真正の逆行列P-1との間に差分dP-1が存在したものとする。このことを式で表すと下記の式となる。
【数４８】
2015113121.tif 000051
ここで、P*dP-1=dEであることに着目すれば、ガウスの消去法を用いることで（マシン精度と同一精度で、しかもはるかに高速で）差分dP-1を求めることができる。その場合には、真正の逆行列P-1は実際に算出されることはなく、差分dP-1は下式により表される。
【数４９】
2015113121.tif 000052
続いて、
【数５０】
2015113121.tif 000053
と置き、そして以上の手順を反復実行して、所与の判定基準が満たされた時点で反復実行を終了する。所与の判定基準としては、例えば、
【数５１】
2015113121.tif 000054
を用いることができる。尚、ここで下付文字のn-th iterationは、n回目の反復実行であることを表している。
【００７６】
本発明によれば、任意の特定の時点における自己汚染の程度を一度の計算ステップで直接的に算出することができ、時間のかかる無用な中間ステップを反復実行することを必要としない。更に、算出値として得られる解は、計算精度の範囲内での厳密性を有する。本発明に係る方法は、宇宙航行機に関する計算及びシミュレーションのためのソフトウェア・ツール・スイートに、その構成要素の１つとして組込むのに適したものである。
【図１】
2015113121.tif 000059
【図２】
2015113121.tif 000060
【図３】
2015113121.tif 000061
【図４】
2015113121.tif 000062