(57)【特許請求の範囲】
【請求項１】
軸対称タンクを球座標で表し、該軸対称タンクの内部の液体を非粘性と仮定して非粘性流体の表面張力を含むスロッシングの基本次モードに関する振動方程式を導き、前記軸対称タンク壁面での前記液体の摩擦力と表面張力との釣り合いを履歴を考慮してモデル化し、前記摩擦力を前記液体の前記軸対称タンクに対する接触角の変化で表す関係式を導き、前記摩擦力の仮想仕事を評価し、この仮想仕事を前記振動方程式に付加した非線形減衰方程式を導き、前記非線形減衰方程式の等価線形方程式に対する減衰比の表現を導くことにより得られた前記減衰比の表現を含む所定のプログラムを格納した記憶手段と、
前記軸対称タンクの寸法、前記液体の密度、前記接触角、前記摩擦力による前記接触角の変化、及び前記液体の前記軸対称タンクの壁面における振幅の少なくとも１つの値を入力する入力手段と、
前記記憶手段に格納された前記プログラムを読み込んで実行し、前記入力手段に入力された前記少なくとも１つの値に基づいて前記減衰比の表現の値を計算し、この値が閾値を超えたかどうかを判定する演算手段と、
前記演算手段で計算した前記減衰比の値と前記判定の結果を出力する出力手段と
を含むことを特徴とするスロッシング減衰比の計算装置。
【請求項２】
前記軸対称タンクは球形であり、前記軸対称タンクの寸法は、該軸対称タンクの半径であることを特徴とする請求項１記載のスロッシング減衰比の計算装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、宇宙機のタンク等の低重力場に配置されるタンク内の貯蔵液体におけるスロッシング振幅の減衰を貯蔵液体のタンクに対する接触角履歴を考慮して予測するスロッシング減衰比の予測方法及び該スロッシング減衰比の計算装置に関する。
【背景技術】
【０００２】
一般に、液体を貯蔵するタンクについてタンク構造設計する場合は、タンク内で生じる液体(貯蔵液体)のスロッシングに伴うスロッシング動液圧による荷重の正確な予測が必要不可欠である。例えば、宇宙衛星等の宇宙機に搭載される推進薬タンクのような液体貯蔵用のタンクでは、該タンク内で液体のスロッシングが生じると、宇宙機自体の姿勢制御に影響を及ぼす虞がある。そのために、この種の宇宙機に搭載される液体貯蔵用のタンクにおいても、貯蔵された液体のスロッシングについての減衰の予測が求められる。
【０００３】
なお、地表重力場に設けられたタンクにおいては、タンクに貯蔵された液体に働く表面張力は、該液体に作用する重力に比して大幅に小さいために、該タンク内の液体の液面はほぼフラットで、その外周縁部にわずかに表面張力の影響によるメニスカスが形成されるに過ぎない。よって、地表重力場におけるタンク内液体のスロッシングの減衰比を予測する場合は、該液体の表面張力はほぼ無視することができる。
【０００４】
しかし、宇宙空間は低重力場であるため、宇宙機のタンク内に貯蔵された液体については、重力の影響が小さくなることに伴って表面張力の影響が大となり、形成されるメニスカスも大きくなる。したがって
気体と液体との界面（気液界面）
固体と液体との界面（固液界面）
固体と気体との界面（固気界面）
での界面張力によって決まる、液面とタンク壁面との接触角が重要になる。このため減衰は、粘性境界層だけでなく、接触角履歴によっても生じることになる。
【０００５】
ここで、接触角履歴とは、液面とタンク壁面とのなす接触角が、摩擦により、液面上昇時には静的平衡時の値よりも大きく、液面下降時には小さくなる現象である。液面が同じ位置でも上昇中か下降中かによって接触角が異なり、接触角が履歴に依存するので、接触角履歴と称される。
【０００６】
従来、接触角履歴のモデル化方法として、液面を張力膜と考え膜の周辺にばね支持境界条件を課す方法が慣用されている。非特許文献１では、式（１６）とその下の説明のように、液面変位の半径方向勾配が変位に比例する形で接触角履歴を表し、特許文献１に挙げられた文献４，９でも同様な式で表されていることを記している。
【０００７】
しかしながら、このようなモデル化方法では、摩擦によるエネルギ散逸を考慮できないため、接触角履歴による減衰比を予測することはできなかった。非特許文献２では、ｐ．１７１，式（５６ａ）とその下の説明のように、履歴（ヒステリシス）による減衰、エネルギ散逸が説明できないことを述べている。
【０００８】
一方、非特許文献３には、低重力における軸対称容器内で軸方向に加振される液体の基本次モードのスロッシングに関する非減衰の振動方程式が開示されている。
【先行技術文献】
【非特許文献】
【０００９】
【非特許文献１】Peterson, L.D., Crawley, E.F., and Hansman, R.J., “宇宙機の力学特性と連成した非線形の燃料スロッシング（Nonlinear Fluid Slosh Coupled to the Dynamics of a Spacecraft）”，エイアイエイエイ（アメリカン  インスティチュート  オブ  エアロノーティックス  アンド  アストロノーティックス；米国航空宇宙学会）  ジャ−ナル（AIAA Journal）(宇宙航行機とロケットの雑誌、論文集) l, Vol.27, No. 9, 1989, pp. 1230-1240
【非特許文献２】Ibrahim, R. A., Pilipchuk, V. N. and Ikeda, T., 2001, “液体スロッシングの力学特性に関する最近の進展（Recent Advances in Liquid Sloshing Dynamics）”，アプライド  メカニクス  レビュウ（Applied Mechanics Reviews）（応用力学のレビュウ), Vol.54, 2001, pp. 133-199
【非特許文献３】内海（Utsumi, M.）, “軸方向に加振される軸対称容器内の低重力スロッシング（Low-gravity Sloshing in an Axisymmetrical Container Excited in the Axial Direction）”， アスメ（アメリカン  ソサイエティ  オブ  メカニカル  エンジニアズ；米国機械学会)ジャ−ナル  オブ  アプライド  メカニクス（ASME Journal of Applied Mechanics）(応用力学の雑誌、論文集)”， Vol. 67, June 2000, pp. 344-354
【発明の概要】
【発明が解決しようとする課題】
【００１０】
従来の接触角履歴によるスロッシングのモデル化方法においては、摩擦によるエネルギ散逸を考慮できないため、接触角履歴による減衰比を予測することはできなかった。このため、エネルギ散逸も考慮した接触角履歴によるスロッシング減衰比の予測方法及び該スロッシング減衰比の計算装置が求められていた。
【課題を解決するための手段】
【００１１】
本願に係るスロッシング減衰比の予測方法の発明は、軸対称タンクを球座標で表し、該一般軸対称タンクの内部の液体を非粘性と仮定して非粘性流体の表面張力を含むスロッシングの基本次モードに関する振動方程式を導き、前記軸対称タンク壁面での前記液体の摩擦力と表面張力との釣り合いを履歴を考慮してモデル化し、前記摩擦力を前記液体の前記軸対称タンクに対する接触角の変化で表す関係式を導き、前記摩擦力の仮想仕事を評価し、この仮想仕事を前記振動方程式に付加した非線形減衰方程式を導き、前記非線形減衰方程式の等価線形方程式に対する減衰比の表現を導くものである。
【００１２】
前記軸対称タンクは球形であり、前記軸対称タンクの半径、前記液体の密度、前記接触角、前記摩擦力による前記接触角の変化、及び前記液体の前記軸対称タンクの壁面における振幅の少なくとも１つの値に基づいて前記減衰比の値を算定することが好ましい。
【００１３】
本願に係る計算装置の発明は、軸対称タンクを球座標で表し、該一般軸対称タンクの内部の液体を非粘性と仮定して非粘性流体の表面張力を含むスロッシングの基本次モードに関する振動方程式を導き、前記軸対称タンク壁面での前記液体の摩擦力と表面張力との釣り合いを履歴を考慮してモデル化し、前記摩擦力を前記液体の前記軸対称タンクに対する接触角の変化で表す関係式を導き、前記摩擦力の仮想仕事を評価し、この仮想仕事を前記振動方程式に付加した非線形減衰方程式を導き、前記非線形減衰方程式の等価線形方程式に対する減衰比の表現を導くことにより得られた前記減衰比の表現を含む所定のプログラムを格納した記憶手段と、前記軸対称タンクの寸法、前記液体の密度、前記接触角、前記摩擦力による前記接触角の変化、及び前記液体の前記軸対称タンクの壁面における振幅の少なくとも１つの値を入力する入力手段と、前記記憶手段に格納された前記プログラムを読み込んで実行し、前記入力手段に入力された前記少なくとも１つの値に基づいて前記減衰比の表現の値を計算し、この値が前記閾値を超えたかどうかを判定する演算手段と、前記演算手段で計算した前記減衰比の値と前記判定の結果を出力する出力手段とを含むものである。
【００１４】
前記軸対称タンクは球形であり、前記軸対称タンクの寸法は、該軸対称タンクの半径であることが好ましい。
【発明の効果】
【００１５】
本発明によると、エネルギ散逸も考慮した接触角履歴によるスロッシング減衰比を予測することができる。
【図面の簡単な説明】
【００１６】
【図１】軸対称タンクを計算モデルとする場合の座標系を示す図である。
【図２】タンク壁面における界面張力のつりあいを示し、ａ）静的平衡時、ｂ）タンク壁面での液面速度が上向きの場合、ｃ）タンク壁面での液面速度が下向きの場合を示す図である。
【図３】等価線形減衰比ζｅｑを示す図である。
【図４】非線形摩擦減衰項の係数Ｃｆを示す図である。
【図５】質量パラメータＭｓを示す図である。
【図６】Ｃｆ／Ｍｓを示す図である。
【図７】固有振動数ω／２πを示す図である。
【図８】メニスカスの面積を示す図である。
【図９】低いボンド数でのメニスカスの接触線と中心のz座標を示す図である。
【図１０】粘性境界層による減衰比を示す図である。
【図１１】接触角履歴によるスロッシング減衰比を予測する計算装置の構成を示す図である。
【図１２】計算装置における一連の動作の流れを示すフローチャートである。
【発明を実施するための形態】
【００１７】
以下、本発明に係る接触角履歴によるスロッシング減衰の予測方法及び該スロッシング減衰比の計算装置の実施の形態について、図面を参照して詳細に説明する。
【００１８】
〔１．実施の形態の構成〕
本実施の形態は、以下のような順序で構成されている。最初に、非特許文献３の方法で、基本次モードに関する非減衰の振動方程式を導き（２．２節）、この方程式に接触角履歴による減衰項を次の手順で導入する。
（１）摩擦力と界面張力のつりあいを、履歴を考慮して記述する正確なモデル化を行い、摩擦力を接触角変化で表す関係式を導く（２．３節）。
（２）摩擦力（非保存力）の仮想仕事を評価し、この仮想仕事を非減衰方程式に付加する（２．４節）。
（３）履歴非線形により非線形減衰となるので、等価線形系の減衰比（等価減衰比）を算定する（２．５節）。
【００１９】
実施例（３節）の結果、重力が低下すると等価減衰比は、履歴による非線形によって、粘性減衰比よりも顕著に増大すること、このために、本明細書で求めた接触角履歴による減衰比を省略すると実際の減衰比（粘性と接触角履歴による減衰の和）を１／３以下にまで過小評価し得ることが分った。
【００２０】
さらに、上述のような非線形減衰を考慮したスロッシング減衰比について、計算を実行して予測を可能にする計算装置を具体例として示す（４節）。また、付録Ａとして式（３２）の導出を示す。
【００２１】
〔２．計算方法〕
〔２．１  計算モデル〕
人工衛星、宇宙航行機に多用される任意軸対称タンク（円筒タンク以外を称し、球形タンクが例として頻用される）を対象とし、図１のような軸対称タンク１内スロッシングの計算モデルを考える。
Ｖ：液体領域
Ｗ：タンク壁面
Ｍ：静的平衡時の液面（メニスカス）
Ｆ：振動している液面
ζ：液面の振動変位
【００２２】
メニスカスは、低重力場では表面張力により強く湾曲した軸対称面となり、表面張力の無視できる地表重力場ではタンク対称軸（ｚ軸）に垂直な平面となる。液体２の運動は非圧縮完全流体の渦なし流れ、タンク１は剛体と仮定する。液面の振動振幅は小さいとして液面境界条件に線形理論を用いる。すなわち、接触角履歴による非線形減衰以外には、線形理論を用いる。
【００２３】
〔２．２  非減衰のモード方程式〕
本明細書では、接触角履歴を引き起こす摩擦力による減衰項を、摩擦力のなす仮想仕事から導く。このため、まず、減衰を導入する前の非減衰モード方程式を、エネルギ原理に基づく変分形で導いておく必要がある。接触角条件を決める界面エネルギを考慮したハミルトンの原理は、次式によって与えられる。
【数１】
0005871121.tif 000002
【００２４】
［  ］の中がラグランジュアンを表す。液圧の項以外、すなわち、気圧と界面エネルギによるポテンシャルエネルギの項を省略すると、液圧がラグランジュアン密度に等しいことに基づく通常の（低重力でない）スロッシング問題のハミルトン原理に帰着する。
【００２５】
ハミルトンの原理［式（１）］を、解析的方法によりモード離散化するため、図１のような球座標系
【数２】
0005871121.tif 000003
を導入し、メニスカスＭ，振動液面Ｆ，タンク壁面Ｗの動径座標を角座標の関数として次のように表す。
【数３】
0005871121.tif 000004
【００２６】
球座標の原点Ｏはメニスカスとタンク壁面との接触交線でタンク壁面に接する円錐の頂点である。このため、
ζ：Ｒ方向に設定した液面変位
は適合条件（タンク壁面で液面変位が壁面を貫いたり壁面から離れたりせず、壁面に接する条件）を満たす。これは、本解析法独自の工夫で、本発明の、接触角履歴のモデル化にも便利である。
【００２７】
タンク壁面のｒ座標のｚ微分がメニスカスの接触交線で正の場合には、球座標の原点は、図１と対照的にタンクの下側になる。従って、球座標と円筒座標の関係は
【数４】
0005871121.tif 000005
【００２８】
減衰比の予測には、タンクの加振加速度のない自由振動を考えれば十分であり、加振加速度項を省略して、液圧ｐｌを圧力方程式から次のように表す。
【数５】
0005871121.tif 000006
【００２９】
式（６）を式（１）に代入し、液体領域Ｖが液面変位を介して変動することに注意して変分計算を実行すると、次式が導かれる。
【数６】
0005871121.tif 000007
【００３０】
式（７）より、変分
【数７】
0005871121.tif 000008
の任意独立性より下記のような支配方程式系が導かれる：
【数８】
0005871121.tif 000009
【００３１】
変分原理（７）を球座標で表し、モード離散化して、変分形モード方程式を次の形に得る。
【数９】
0005871121.tif 000010
【００３２】
Ｍｓ、Ｋｓはスロッシングの質量、剛性パラメータで、固有振動数は次式によって定まる。
【数１０】
0005871121.tif 000011
【００３３】
モード方程式の解であるモード座標ｑを用いて、速度ポテンシャル、液面変位は次式によって表される。
【数１１】
0005871121.tif 000012
【００３４】
任意軸対称タンクについて特性関数を解析的に導出可能とし、計算を大幅に効率化したことが非特許文献３の特長である。
【００３５】
〔２．３  摩擦力を接触角の静的接触角からの変化で表す関係〕
界面張力の静的つりあい式は、図２（ａ）より、次式によって与えられる。
【数１２】
0005871121.tif 000013
【００３６】
式（１７）は下記を意味するため、静的な接触角条件と呼ばれる。
【数１３】
0005871121.tif 000014
【００３７】
スロッシング時の動的な接触角条件は、摩擦による接触角履歴がないとき、上で導入した変分原理（１）の停留条件Ｅ５＝０［式（１２）参照］である。これらの静的および動的な接触角条件の比較によって、摩擦による接触角履歴がないとき、次の関係が成立することが分かる。
【数１４】
0005871121.tif 000015
【００３８】
スロッシング時の界面張力と摩擦力の釣り合い式は、図２（ｂ），（ｃ）双方の場合をまとめて次のように表される。
【数１５】
0005871121.tif 000016
【００３９】
宇宙機の液体推進薬スロッシングでは、静的接触角は小さく（５度前後）、その変化も小さいので、式（１８）は次のように近似できる。
【数１６】
0005871121.tif 000017
【００４０】
式（１９）においてｓｇｎを含まない項は、静止時の釣り合い式（１７）より消える。従って、摩擦力を接触角変化で表す次の関係式が導かれる。
【数１７】
0005871121.tif 000018
【００４１】
〔２．４  接触角履歴による減衰項の導出〕
摩擦力による仮想仕事δＷｆは、
【数１８】
0005871121.tif 000019
を接触線に沿って積分することにより、次式によって計算される：
【数１９】
0005871121.tif 000020
【００４２】
式（２０）を代入して
【数２０】
0005871121.tif 000021
【００４３】
式（２２）のｓｇｎ関数の部分は、液面変位のモード展開式（１６）を代入すると
【数２１】
0005871121.tif 000022
【００４４】
式（２３）の下線部は、上向き液面変位のモード関数の、タンク壁面上の位置
【数２２】
0005871121.tif 000023
での値である。この値が１となるように定数ｃｋを規格化し、式（２３）を次のように変形する。
【数２３】
0005871121.tif 000024
【００４５】
式（２４）を式（２２）に代入すると下記のようになる。
【数２４】
0005871121.tif 000025
【００４６】
式（１６）を代入して
【数２５】
0005871121.tif 000026
【００４７】
積分は、
【数２６】
0005871121.tif 000027
となり、次式を得る。
【数２７】
0005871121.tif 000028
【００４８】
式（２６）を式（１４）の左辺に加算することによって、接触角履歴による非線形減衰挙動を支配する振動方程式が次のように導かれる。
【数２８】
0005871121.tif 000029
【００４９】
式（２８）は次のように書ける。
【数２９】
0005871121.tif 000030
ここで
【数３０】
0005871121.tif 000031
【００５０】
〔２．５  等価線形減衰比〕
式（２９）に関して、次のような等価線形化方程式
【数３１】
0005871121.tif 000032
を考え、付録Ａに示す方法で、次のように等価線形減衰比を求める。
【数３２】
0005871121.tif 000033
【００５１】
式（２７），（３２）より、等価線形減衰比は、下記に比例する：
【数３３】
0005871121.tif 000034
【００５２】
〔３．実施例〕
〔３．１  減衰比の評価結果〕
下記のパラメータに関して数値計算に用いた。
【数３４】
0005871121.tif 000035
【００５３】
重力加速度は、ボンド数
【数３５】
0005871121.tif 000036
によって定めた。これは、重力と表面張力の効果の度合を表わす無次元数である。
【００５４】
図３に、等価線形減衰比を示す。横軸は充填率（液体の体積のタンクの体積に対する比）である。図３より、下記の２点が分る：
（１）  等価線形減衰比は、ボンド数が低下すると増加する。
（２）  等価線形減衰比は、特に小さいボンド数１，０．１に関し、充填率の増加と共に単調増加するようになる。
【００５５】
これらの結果（１），（２）の理由を調べるため、式（３２）右辺に現れる諸パラメータの変化を次のように図４，５，６，７に示す。
図４：非線形摩擦減衰項の係数Ｃｆ
図５：質量パラメータＭｓ
図６：Ｃｆ／Ｍｓ
図７：固有振動数ω／２π
【００５６】
〔３．２  結果（１）に関する考察〕
まず、結果（１）の理由について考察する。図４，５より、ボンド数Ｂｏの低下に伴う下記の傾向が観察される。
（ａ）概ね０．５より高い充填率では、ＣｆはＭｓと異なり単調に減少する。
（ｂ）概ね０．５より低い充填率では、ＭｓはＣｆと異なり単調に増加する。
【００５７】
その結果、図６のように、Ｃｆ／Ｍｓはボンド数の低下に伴い減少する（ボンド数１００以上の高い充填率でのみ例外的傾向が見られる）。しかし、式（３２）分母に表れる固有振動数ωは、図７のように、ボンド数の減少に伴いＣｆ／Ｍｓよりもさらに著しく低下する。その結果、等価線形減衰比は、図３のように、ボンド数が減少すると増加する。このように、３．１節の結果（１）がボンド数の減少による固有振動数の低下に起因することが分った。
【００５８】
次の２つの節では、図４，５に示されたＣｆ、Ｍｓの変化について議論しておく。
【００５９】
〔３．３  Ｃｆの変化について議論〕
図４に示されたＣｆの変化は、式（２７）より
【数３６】
0005871121.tif 000037
の変化を表している。これは、式（２）、図１より、メニスカスの接触線の半径である。メニスカスの接触線の半径は、ボンド数無限大、すなわち表面張力がなくメニスカスが平面のとき、充填率０．５で最大でタンク半径に等しい。充填率一定でボンド数が低下すると、メニスカスの下に凸の曲がりが顕著になるために、メニスカスの接触線は上昇する。従って、図４のように、充填率０．５以上のときは、ボンド数が低下すると、接触線がタンクの頂上に近づき接触線の半径は減少する。また、メニスカスの接触線の半径を最大にする充填率（接触線のｚ座標がタンク半径に等しくなる充填率）は、ボンド数の低下とともに減少する。
【００６０】
〔３．４  Ｍｓの変化について議論〕
図５に示された質量パラメータＭｓの変化は、図８に示されたメニスカスの面積に似た変化を呈する。この理由は次のように説明される。質量パラメータＭｓは運動エネルギを表し、運動エネルギはグリーンの定理とタンク壁面上の境界条件
【数３７】
0005871121.tif 000038
を用いてメニスカス上の面積分：
【数３８】
0005871121.tif 000039
で表される。この積分値はメニスカスの面積に大きく依存する。従って、質量パラメータＭｓとメニスカス面積のボンド数、充填率に対する依存性は類似してくる。
【００６１】
〔３．５  結果（２）に関する考察〕
３．１節の結果（２）の理由について考察する。図７より、ボンド数が小さい１．０，０．１のとき、固有振動数ωの充填率依存性は弱い。従って、式（３２）によって与えられる等価線形減衰比は、Ｃｆ／Ｍｓと同様な充填率依存性を有する（図３，６参照）。従って、結果（２）の理由は、Ｃｆ／Ｍｓの変動の原因を考察することによって調べられる。
【００６２】
図６より、Ｃｆ／Ｍｓは低い充填率よりも高い充填率で大きく、図４より、ボンド数が小さい１．０，０．１のとき、Ｃｆは高い充填率よりも低い充填率で大きい。従って、３．１節の結果（２）は、ボンド数が小さい１．０，０．１のときＭｓが高い充填率より低い充填率で大きくなることに基づく（図５参照）。３．４節で記したように、質量パラメータＭｓはメニスカスの面積と強い相関をもつ。従って、３．１節の結果（２）の理由は、メニスカスの面積が高い充填率より低い充填率で大きくなることに基づく（図８参照）。
【００６３】
メニスカスの面積の指標は、メニスカスの接触線と中心のｚ座標の差である。図９に、これらのｚ座標をボンド数１，０．１の場合について示した。ボンド数が１から０．１まで低下する際の低い充填率でのメニスカスの面積の増加は、主として、メニスカス中心の下降ではなく、接触線の上昇に起因することが分かる。
【００６４】
〔３．６  粘性減衰比との比較〕
図１０に、粘性減衰比ζｖｉｓを示す。粘性係数はμ＝０．００１１Ｎｓ／ｍ２である。図３と図１０を比較することによって、ボンド数の低下と共に増加する傾向が、接触角履歴による等価減衰比の方が粘性減衰比よりも強いことが確認できる。この理由は次のように考察できる。粘性による減衰比は、固有振動数の逆数に比例する。接触角履歴による等価線形減衰比は、固有振動数の逆数ではなく、その２乗に比例する。固有振動数は、図７のように、ボンド数が低下すると減少する。このため、接触角履歴による減衰比は、粘性減衰比に比べて、ボンド数の低下と共に増加する傾向が強くなる。
【００６５】
接触角履歴による減衰比と粘性減衰比の固有振動数に対する依存性の違いは、履歴非線形に基づく（付録の最後の段落参照）。
【００６６】
ボンド数の低下に伴い接触角履歴による減衰比の方が粘性減衰比よりも顕著に増加することは、ボンド数低下と共に、接触角履歴による減衰比の
実減衰比＝接触角履歴による減衰比ζｅｑ＋粘性減衰比ζｖｉｓ
に対する寄与が増大することを意味する。例えば、ボンド数１のとき、
【数３９】
0005871121.tif 000040
は液体充填率０．６，０．９でそれぞれ２．３，３．３に達する。このように、実際のトータル減衰比は、本発明で評価した接触角履歴による減衰比を考慮しなければ、かなり過小評価される。
【００６７】
〔４．計算装置〕
上述した接触角履歴を考慮したスロッシング減衰比の予測は、図１１に示すような計算装置１０によって実現することができる。この計算装置１０は、ＣＰＵ、ＤＳＰの如き演算部１１、ＲＡＭ、ＲＯＭ、ハードディスクの如き記憶部１２、ＬＣＤ、プリンタの如き出力部１３、キーボード、マウスの如き入力部１４を含み、例えばパーソナルコンピュータを利用することができる。
【００６８】
図１２に示す計算装置１１の一連の動作は、記憶部１２に格納されたスロッシング減衰比算定プログラム１２ａを演算部１１が読み出して実行することにより実現される。このスロッシング減衰比算定プログラム１２ａは、前述のような手順によって得られたスロッシング減衰比の表現を含んでいる。
【００６９】
最初のステップＳ１においては、モデルを設定する。ここでは、球形の軸対称タンクを想定し、前記軸対称タンクの半径、前記液体の密度、前記接触角、前記摩擦力による前記接触角の変化、及び前記液体の前記軸対称タンクの壁面における振幅の少なくとも１つの値によりモデルを設定するものとする。入力部１４は、これらの少なくとも１つの値を入力値として受け取る。演算部１１は、入力部１４が受け取った入力値を記憶部１２に格納する。
【００７０】
ステップＳ２においては、演算部１１は、記憶部１２に格納された入力値を読み出し、その数値計算部１１ａにおいて、この入力値に基づいてスロッシング減衰比の表現を用いてこのスロッシング減衰比の値を数値計算する。演算部１１は、得られたスロッシング減衰比の値を記憶部１２に格納する。
【００７１】
ステップＳ３においては、演算部１１は、記憶部１２に格納されたスロッシング減衰比の値と、同じく記憶部１２に格納された所定の閾値１２ｂとを読み出す。演算部１１は、その判定部１１ｂにおいて、スロッシング減衰比の値が閾値１２ｂを超えない場合にはＯＫと判定して一連のステップを終了する。一方、スロッシング減衰比の値が閾値を越えた場合にはＮＧとして判定して前のステップＳ１のモデル設定に手順を戻す。なお、閾値１２ｂは、入力部１４を介して設定することができる。
【００７２】
このような一連の工程において、スロッシング減衰比が所定の閾値内に収まるまでモデル設定、数値計算、判定のループを繰り返すことにより、閾値内に収まるモデル設定を可能としている。また、前述のスロッシング減衰比の表現を利用することにより、エネルギ散逸を考慮した精度の高いモデル設定を可能としている。
【００７３】
なお、このようなスロッシング減衰比の算定は、記憶部１２に格納したスロッシング減衰比算定プログラム１２ａのような、前述のスロッシング減衰比の表現を含み、モデル設定、数値計算、判定のステップを有するプログラムによっても提供することができる。
【００７４】
〔付録Ａ  式（３２）の導出〕
非線形方程式（２９）の解を次のようにおく。
【数４０】
0005871121.tif 000041
ｑの時間微分が
【数４１】
0005871121.tif 000042
となる条件を課すことより、次式が導かれる。
【数４２】
0005871121.tif 000043
【００７５】
式（Ａ２）を時間で微分すると、次のようになる。
【数４３】
0005871121.tif 000044
【００７６】
式（Ａ１），  （Ａ２），  （Ａ４）  を式（２９）に代入すると次式が得られる。
【数４４】
0005871121.tif 000045
【００７７】
式（Ａ３），  （Ａ５）をまとめて記すと
【数４５】
0005871121.tif 000046
ここで、次のようにおいた。
【数４６】
0005871121.tif 000047
【００７８】
式（Ａ６）を
【数４７】
0005871121.tif 000048
について解くと
【数４８】
0005871121.tif 000049
【００７９】
Ａ、Φは徐々に変動する時間関数であるため、式（Ａ７）の右辺は１周期の時間平均として近似できる：
【数４９】
0005871121.tif 000050
【００８０】
式（Ａ８）の第２式の被積分関数はΨについて奇関数であることに留意して積分を実行すると
【数５０】
0005871121.tif 000051
が得られる。
【００８１】
式（３１）の特性方程式
【数５１】
0005871121.tif 000052
の根の実部
【数５２】
0005871121.tif 000053
より、等価線形系の減衰振動の振幅は
【数５３】
0005871121.tif 000054
となり、次の関係を満たす。
【数５４】
0005871121.tif 000055
【００８２】
式（Ａ９），（Ａ１２）から
【数５５】
0005871121.tif 000056
を消去すれば式（３２）が導かれる。すなわち、振幅変動が式（Ａ９）によって与えられる非線形系の等価線形減衰比は、式（３２）によって与えられる。
【００８３】
減衰比の固有振動数ωへの依存性を、線形減衰の場合と比較する。線形の場合の依存性を調べるため、次の変更：
【数５６】
0005871121.tif 000057
を行うと、次式が得られる
【数５７】
0005871121.tif 000058
【００８４】
式（Ａ１３）を
【数５８】
0005871121.tif 000059
と比較すると、線形減衰比は次のように固有振動数の逆数に比例することが分かる。
【数５９】
0005871121.tif 000060
【００８５】
接触角履歴を有する非線形系の等価線形減衰比は、式（３２）のように、固有振動数ではなく、その２乗の逆数に比例し、線形減衰比よりも固有振動数依存性は強まることが分かる。
【００８６】
なお、上述の実施の形態は、本発明の一具体例を示すものであり、本発明を限定するものではない。
【符号の説明】
【００８７】
１０  計算装置
１１  演算部
１２  記憶部
１３  出力部
１４  入力部
【図１】
0005871121.tif 000061
【図２】
0005871121.tif 000062
【図３】
0005871121.tif 000063
【図４】
0005871121.tif 000064
【図５】
0005871121.tif 000065
【図６】
0005871121.tif 000066
【図７】
0005871121.tif 000067
【図８】
0005871121.tif 000068
【図９】
0005871121.tif 000069
【図１０】
0005871121.tif 000070
【図１１】
0005871121.tif 000071
【図１２】
0005871121.tif 000072