"軸対称タンクを球座標で表し、該軸対称タンクの内部の液体を非粘性と仮定して非粘性流体の表面張力を含むスロッシングの基本次モードに関する振動方程式を導き、前記軸対称タンク壁面での前記液体の摩擦力と表面張力との釣り合いを履歴を考慮してモデル化し、前記摩擦力を前記液体の前記軸対称タンクに対する接触角の変化で表す関係式を導き、前記摩擦力の仮想仕事を評価し、この仮想仕事を前記振動方程式に付加した非線形減衰方程式を導き、前記非線形減衰方程式の等価線形方程式に対する減衰比の表現を導くことを特徴とするスロッシング減衰比の予測方法。前記軸対称タンクは球形であり、前記軸対称タンクの半径、前記液体の密度、前記接触角、前記摩擦力による前記接触角の変化、及び前記液体の前記軸対称タンクの壁面における振幅の少なくとも1つの値に基づいて前記減衰比の値を算定することを特徴とする請求項1記載のスロッシング減衰比の予測方法。軸対称タンクを球座標で表し、該軸対称タンクの内部の液体を非粘性と仮定して非粘性流体の表面張力を含むスロッシングの基本次モードに関する振動方程式を導き、前記軸対称タンク壁面での前記液体の摩擦力と表面張力との釣り合いを履歴を考慮してモデル化し、前記摩擦力を前記液体の前記軸対称タンクに対する接触角の変化で表す関係式を導き、前記摩擦力の仮想仕事を評価し、この仮想仕事を前記振動方程式に付加した非線形減衰方程式を導き、前記非線形減衰方程式の等価線形方程式に対する減衰比の表現を導くことにより得られた前記減衰比の表現を含む所定のプログラムを格納した記憶手段と、前記軸対称タンクの寸法、前記液体の密度、前記接触角、前記摩擦力による前記接触角の変化、及び前記液体の前記軸対称タンクの壁面における振幅の少なくとも1つの値を入力する入力手段と、前記記憶手段に格納された前記プログラムを読み込んで実行し、前記入力手段に入力された前記少なくとも1つの値に基づいて前記減衰比の表現の値を計算し、この値が前記閾値を超えたかどうかを判定する演算手段と、前記演算手段で計算した前記減衰比の値と前記判定の結果を出力する出力手段とを含むことを特徴とするスロッシング減衰比の計算装置。前記軸対称タンクは球形であり、前記軸対称タンクの寸法は、該軸対称タンクの半径であることを特徴とする請求項3記載のスロッシング減衰比の計算装置。従来の接触角履歴によるスロッシングのモデル化方法においては、摩擦によるエネルギ散逸を考慮できないため、接触角履歴による減衰比を予測することはできなかった。このため、エネルギ散逸も考慮した接触角履歴によるスロッシング減衰比の予測方法及び該スロッシング減衰比の計算装置が求められていた。本願に係るスロッシング減衰比の予測方法の発明は、軸対称タンクを球座標で表し、該一般軸対称タンクの内部の液体を非粘性と仮定して非粘性流体の表面張力を含むスロッシングの基本次モードに関する振動方程式を導き、前記軸対称タンク壁面での前記液体の摩擦力と表面張力との釣り合いを履歴を考慮してモデル化し、前記摩擦力を前記液体の前記軸対称タンクに対する接触角の変化で表す関係式を導き、前記摩擦力の仮想仕事を評価し、この仮想仕事を前記振動方程式に付加した非線形減衰方程式を導き、前記非線形減衰方程式の等価線形方程式に対する減衰比の表現を導くものである。前記軸対称タンクは球形であり、前記軸対称タンクの半径、前記液体の密度、前記接触角、前記摩擦力による前記接触角の変化、及び前記液体の前記軸対称タンクの壁面における振幅の少なくとも1つの値に基づいて前記減衰比の値を算定することが好ましい。本願に係る計算装置の発明は、軸対称タンクを球座標で表し、該一般軸対称タンクの内部の液体を非粘性と仮定して非粘性流体の表面張力を含むスロッシングの基本次モードに関する振動方程式を導き、前記軸対称タンク壁面での前記液体の摩擦力と表面張力との釣り合いを履歴を考慮してモデル化し、前記摩擦力を前記液体の前記軸対称タンクに対する接触角の変化で表す関係式を導き、前記摩擦力の仮想仕事を評価し、この仮想仕事を前記振動方程式に付加した非線形減衰方程式を導き、前記非線形減衰方程式の等価線形方程式に対する減衰比の表現を導くことにより得られた前記減衰比の表現を含む所定のプログラムを格納した記憶手段と、前記軸対称タンクの寸法、前記液体の密度、前記接触角、前記摩擦力による前記接触角の変化、及び前記液体の前記軸対称タンクの壁面における振幅の少なくとも1つの値を入力する入力手段と、前記記憶手段に格納された前記プログラムを読み込んで実行し、前記入力手段に入力された前記少なくとも1つの値に基づいて前記減衰比の表現の値を計算し、この値が前記閾値を超えたかどうかを判定する演算手段と、前記演算手段で計算した前記減衰比の値と前記判定の結果を出力する出力手段とを含むものである。前記軸対称タンクは球形であり、前記軸対称タンクの寸法は、該軸対称タンクの半径であることが好ましい。本発明によると、エネルギ散逸も考慮した接触角履歴によるスロッシング減衰比を予測することができる。以下、本発明に係る接触角履歴によるスロッシング減衰の予測方法及び該スロッシング減衰比の計算装置の実施の形態について、図面を参照して詳細に説明する。1.実施の形態の構成本実施の形態は、以下のような順序で構成されている。最初に、非特許文献3の方法で、基本次モードに関する非減衰の振動方程式を導き2.2節、この方程式に接触角履歴による減衰項を次の手順で導入する。1摩擦力と界面張力のつりあいを、履歴を考慮して記述する正確なモデル化を行い、摩擦力を接触角変化で表す関係式を導く2.3節。2摩擦力非保存力の仮想仕事を評価し、この仮想仕事を非減衰方程式に付加する2.4節。3履歴非線形により非線形減衰となるので、等価線形系の減衰比等価減衰比を算定する2.5節。実施例3節の結果、重力が低下すると等価減衰比は、履歴による非線形によって、粘性減衰比よりも顕著に増大すること、このために、本明細書で求めた接触角履歴による減衰比を省略すると実際の減衰比粘性と接触角履歴による減衰の和を1/3以下にまで過小評価し得ることが分った。さらに、上述のような非線形減衰を考慮したスロッシング減衰比について、計算を実行して予測を可能にする計算装置を具体例として示す4節。また、付録Aとして式32の導出を示す。2.計算方法2.1計算モデル人工衛星、宇宙航行機に多用される任意軸対称タンク円筒タンク以外を称し、球形タンクが例として頻用されるを対象とし、図1のような軸対称タンク1内スロッシングの計算モデルを考える。V:液体領域W:タンク壁面M:静的平衡時の液面メニスカスF:振動している液面ζ:液面の振動変位メニスカスは、低重力場では表面張力により強く湾曲した軸対称面となり、表面張力の無視できる地表重力場ではタンク対称軸z軸に垂直な平面となる。液体2の運動は非圧縮完全流体の渦なし流れ、タンク1は剛体と仮定する。液面の振動振幅は小さいとして液面境界条件に線形理論を用いる。すなわち、接触角履歴による非線形減衰以外には、線形理論を用いる。2.2非減衰のモード方程式本明細書では、接触角履歴を引き起こす摩擦力による減衰項を、摩擦力のなす仮想仕事から導く。このため、まず、減衰を導入する前の非減衰モード方程式を、エネルギ原理に基づく変分形で導いておく必要がある。接触角条件を決める界面エネルギを考慮したハミルトンの原理は、次式によって与えられる。[]の中がラグランジュアンを表す。液圧の項以外、すなわち、気圧と界面エネルギによるポテンシャルエネルギの項を省略すると、液圧がラグランジュアン密度に等しいことに基づく通常の低重力でないスロッシング問題のハミルトン原理に帰着する。ハミルトンの原理[式1]を、解析的方法によりモード離散化するため、図1のような球座標系を導入し、メニスカスM,振動液面F,タンク壁面Wの動径座標を角座標の関数として次のように表す。球座標の原点Oはメニスカスとタンク壁面との接触交線でタンク壁面に接する円錐の頂点である。このため、ζ:R方向に設定した液面変位は適合条件タンク壁面で液面変位が壁面を貫いたり壁面から離れたりせず、壁面に接する条件を満たす。これは、本解析法独自の工夫で、本発明の、接触角履歴のモデル化にも便利である。タンク壁面のr座標のz微分がメニスカスの接触交線で正の場合には、球座標の原点は、図1と対照的にタンクの下側になる。従って、球座標と円筒座標の関係は減衰比の予測には、タンクの加振加速度のない自由振動を考えれば十分であり、加振加速度項を省略して、液圧plを圧力方程式から次のように表す。式6を式1に代入し、液体領域Vが液面変位を介して変動することに注意して変分計算を実行すると、次式が導かれる。式7より、変分の任意独立性より下記のような支配方程式系が導かれる:変分原理7を球座標で表し、モード離散化して、変分形モード方程式を次の形に得る。Ms、Ksはスロッシングの質量、剛性パラメータで、固有振動数は次式によって定まる。モード方程式の解であるモード座標qを用いて、速度ポテンシャル、液面変位は次式によって表される。任意軸対称タンクについて特性関数を解析的に導出可能とし、計算を大幅に効率化したことが非特許文献3の特長である。2.3摩擦力を接触角の静的接触角からの変化で表す関係界面張力の静的つりあい式は、図2aより、次式によって与えられる。式17は下記を意味するため、静的な接触角条件と呼ばれる。スロッシング時の動的な接触角条件は、摩擦による接触角履歴がないとき、上で導入した変分原理1の停留条件E5=0[式12参照]である。これらの静的および動的な接触角条件の比較によって、摩擦による接触角履歴がないとき、次の関係が成立することが分かる。スロッシング時の界面張力と摩擦力の釣り合い式は、図2b,c双方の場合をまとめて次のように表される。宇宙機の液体推進薬スロッシングでは、静的接触角は小さく5度前後、その変化も小さいので、式18は次のように近似できる。式19においてsgnを含まない項は、静止時の釣り合い式17より消える。従って、摩擦力を接触角変化で表す次の関係式が導かれる。2.4接触角履歴による減衰項の導出摩擦力による仮想仕事δWfは、を接触線に沿って積分することにより、次式によって計算される:式20を代入して式22のsgn関数の部分は、液面変位のモード展開式16を代入すると式23の下線部は、上向き液面変位のモード関数の、タンク壁面上の位置での値である。この値が1となるように定数ckを規格化し、式23を次のように変形する。式24を式22に代入すると下記のようになる。式16を代入して積分は、となり、次式を得る。式26を式14の左辺に加算することによって、接触角履歴による非線形減衰挙動を支配する振動方程式が次のように導かれる。式28は次のように書ける。ここで2.5等価線形減衰比式29に関して、次のような等価線形化方程式を考え、付録Aに示す方法で、次のように等価線形減衰比を求める。式27,32より、等価線形減衰比は、下記に比例する:3.実施例3.1減衰比の評価結果下記のパラメータに関して数値計算に用いた。重力加速度は、ボンド数によって定めた。これは、重力と表面張力の効果の度合を表わす無次元数である。図3に、等価線形減衰比を示す。横軸は充填率液体の体積のタンクの体積に対する比である。図3より、下記の2点が分る:1等価線形減衰比は、ボンド数が低下すると増加する。2等価線形減衰比は、特に小さいボンド数1,0.1に関し、充填率の増加と共に単調増加するようになる。これらの結果1,2の理由を調べるため、式32右辺に現れる諸パラメータの変化を次のように図4,5,6,7に示す。図4:非線形摩擦減衰項の係数Cf図5:質量パラメータMs図6:Cf/Ms図7:固有振動数ω/2π3.2結果1に関する考察まず、結果1の理由について考察する。図4,5より、ボンド数Boの低下に伴う下記の傾向が観察される。a概ね0.5より高い充填率では、CfはMsと異なり単調に減少する。b概ね0.5より低い充填率では、MsはCfと異なり単調に増加する。その結果、図6のように、Cf/Msはボンド数の低下に伴い減少するボンド数100以上の高い充填率でのみ例外的傾向が見られる。しかし、式32分母に表れる固有振動数ωは、図7のように、ボンド数の減少に伴いCf/Msよりもさらに著しく低下する。その結果、等価線形減衰比は、図3のように、ボンド数が減少すると増加する。このように、3.1節の結果1がボンド数の減少による固有振動数の低下に起因することが分った。次の2つの節では、図4,5に示されたCf、Msの変化について議論しておく。3.3Cfの変化について議論図4に示されたCfの変化は、式27よりの変化を表している。これは、式2、図1より、メニスカスの接触線の半径である。メニスカスの接触線の半径は、ボンド数無限大、すなわち表面張力がなくメニスカスが平面のとき、充填率0.5で最大でタンク半径に等しい。充填率一定でボンド数が低下すると、メニスカスの下に凸の曲がりが顕著になるために、メニスカスの接触線は上昇する。従って、図4のように、充填率0.5以上のときは、ボンド数が低下すると、接触線がタンクの頂上に近づき接触線の半径は減少する。また、メニスカスの接触線の半径を最大にする充填率接触線のz座標がタンク半径に等しくなる充填率は、ボンド数の低下とともに減少する。3.4Msの変化について議論図5に示された質量パラメータMsの変化は、図8に示されたメニスカスの面積に似た変化を呈する。この理由は次のように説明される。質量パラメータMsは運動エネルギを表し、運動エネルギはグリーンの定理とタンク壁面上の境界条件を用いてメニスカス上の面積分:で表される。この積分値はメニスカスの面積に大きく依存する。従って、質量パラメータMsとメニスカス面積のボンド数、充填率に対する依存性は類似してくる。3.5結果2に関する考察3.1節の結果2の理由について考察する。図7より、ボンド数が小さい1.0,0.1のとき、固有振動数ωの充填率依存性は弱い。従って、式32によって与えられる等価線形減衰比は、Cf/Msと同様な充填率依存性を有する図3,6参照。従って、結果2の理由は、Cf/Msの変動の原因を考察することによって調べられる。図6より、Cf/Msは低い充填率よりも高い充填率で大きく、図4より、ボンド数が小さい1.0,0.1のとき、Cfは高い充填率よりも低い充填率で大きい。従って、3.1節の結果2は、ボンド数が小さい1.0,0.1のときMsが高い充填率より低い充填率で大きくなることに基づく図5参照。3.4節で記したように、質量パラメータMsはメニスカスの面積と強い相関をもつ。従って、3.1節の結果2の理由は、メニスカスの面積が高い充填率より低い充填率で大きくなることに基づく図8参照。メニスカスの面積の指標は、メニスカスの接触線と中心のz座標の差である。図9に、これらのz座標をボンド数1,0.1の場合について示した。ボンド数が1から0.1まで低下する際の低い充填率でのメニスカスの面積の増加は、主として、メニスカス中心の下降ではなく、接触線の上昇に起因することが分かる。3.6粘性減衰比との比較図10に、粘性減衰比ζvisを示す。粘性係数はμ=0.0011Ns/m2である。図3と図10を比較することによって、ボンド数の低下と共に増加する傾向が、接触角履歴による等価減衰比の方が粘性減衰比よりも強いことが確認できる。この理由は次のように考察できる。粘性による減衰比は、固有振動数の逆数に比例する。接触角履歴による等価線形減衰比は、固有振動数の逆数ではなく、その2乗に比例する。固有振動数は、図7のように、ボンド数が低下すると減少する。このため、接触角履歴による減衰比は、粘性減衰比に比べて、ボンド数の低下と共に増加する傾向が強くなる。接触角履歴による減衰比と粘性減衰比の固有振動数に対する依存性の違いは、履歴非線形に基づく付録の最後の段落参照。ボンド数の低下に伴い接触角履歴による減衰比の方が粘性減衰比よりも顕著に増加することは、ボンド数低下と共に、接触角履歴による減衰比の実減衰比=接触角履歴による減衰比ζeq粘性減衰比ζvisに対する寄与が増大することを意味する。例えば、ボンド数1のとき、は液体充填率0.6,0.9でそれぞれ2.3,3.3に達する。このように、実際のトータル減衰比は、本発明で評価した接触角履歴による減衰比を考慮しなければ、かなり過小評価される。4.計算装置上述した接触角履歴を考慮したスロッシング減衰比の予測は、図11に示すような計算装置10によって実現することができる。この計算装置10は、CPU、DSPの如き演算部11、RAM、ROM、ハードディスクの如き記憶部12、LCD、プリンタの如き出力部13、キーボード、マウスの如き入力部14を含み、例えばパーソナルコンピュータを利用することができる。図12に示す計算装置11の一連の動作は、記憶部12に格納されたスロッシング減衰比算定プログラム12aを演算部11が読み出して実行することにより実現される。このスロッシング減衰比算定プログラム12aは、前述のような手順によって得られたスロッシング減衰比の表現を含んでいる。最初のステップS1においては、モデルを設定する。ここでは、球形の軸対称タンクを想定し、前記軸対称タンクの半径、前記液体の密度、前記接触角、前記摩擦力による前記接触角の変化、及び前記液体の前記軸対称タンクの壁面における振幅の少なくとも1つの値によりモデルを設定するものとする。入力部14は、これらの少なくとも1つの値を入力値として受け取る。演算部11は、入力部14が受け取った入力値を記憶部12に格納する。ステップS2においては、演算部11は、記憶部12に格納された入力値を読み出し、その数値計算部11aにおいて、この入力値に基づいてスロッシング減衰比の表現を用いてこのスロッシング減衰比の値を数値計算する。演算部11は、得られたスロッシング減衰比の値を記憶部12に格納する。ステップS3においては、演算部11は、記憶部12に格納されたスロッシング減衰比の値と、同じく記憶部12に格納された所定の閾値12bとを読み出す。演算部11は、その判定部11bにおいて、スロッシング減衰比の値が閾値12bを超えない場合にはOKと判定して一連のステップを終了する。一方、スロッシング減衰比の値が閾値を越えた場合にはNGとして判定して前のステップS1のモデル設定に手順を戻す。なお、閾値12bは、入力部14を介して設定することができる。このような一連の工程において、スロッシング減衰比が所定の閾値内に収まるまでモデル設定、数値計算、判定のループを繰り返すことにより、閾値内に収まるモデル設定を可能としている。また、前述のスロッシング減衰比の表現を利用することにより、エネルギ散逸を考慮した精度の高いモデル設定を可能としている。なお、このようなスロッシング減衰比の算定は、記憶部12に格納したスロッシング減衰比算定プログラム12aのような、前述のスロッシング減衰比の表現を含み、モデル設定、数値計算、判定のステップを有するプログラムによっても提供することができる。付録A式32の導出非線形方程式29の解を次のようにおく。qの時間微分がとなる条件を課すことより、次式が導かれる。式A2を時間で微分すると、次のようになる。式A1,A2,A4を式29に代入すると次式が得られる。式A3,A5をまとめて記すとここで、次のようにおいた。式A6をについて解くとA、Φは徐々に変動する時間関数であるため、式A7の右辺は1周期の時間平均として近似できる:式A8の第2式の被積分関数はΨについて奇関数であることに留意して積分を実行するとが得られる。式31の特性方程式の根の実部より、等価線形系の減衰振動の振幅はとなり、次の関係を満たす。式A9,A12からを消去すれば式32が導かれる。すなわち、振幅変動が式A9によって与えられる非線形系の等価線形減衰比は、式32によって与えられる。減衰比の固有振動数ωへの依存性を、線形減衰の場合と比較する。線形の場合の依存性を調べるため、次の変更:を行うと、次式が得られる式A13をと比較すると、線形減衰比は次のように固有振動数の逆数に比例することが分かる。接触角履歴を有する非線形系の等価線形減衰比は、式32のように、固有振動数ではなく、その2乗の逆数に比例し、線形減衰比よりも固有振動数依存性は強まることが分かる。なお、上述の実施の形態は、本発明の一具体例を示すものであり、本発明を限定するものではない。"