Activate で Modelica を使ってみる（その 13）- 接触系コンポーネント：ElastoGap と MassWithStopAndFriction
Activate で Modelica を使ってみる（その13）
接触系コンポーネント – ElastoGap と MassWithStopAndFriction
今回と次回で Modelica 標準ライブラリの Mechanics - Translational ライブラリに含まれている接触と摩擦に関わるコンポーネントである ElastoGap、SupportFriction、Brake、MassWithStopAndFriction の４つのコンポーネントについて、摩擦や接触を表現する仕組み、その使い分け及び簡単な例題について紹介します。
この４つのコンポーネントのうち接触を表現できるのはElastoGap と MassWithStopAndFriction で、摩擦を考慮できるのは SupportFriction、Brake、MassWithStopAndFriction の３つです。今回は接触機能を取り上げて、ElastoGap と MassWithStopAndFriction の接触の仕組みと簡単なサンプル事例を紹介します。
ElastoGap
ElastoGap は接触境界をペナルティ法で表現した１次元の接触要素です。ElastoGap で最低限必要となる入力はバネ定数（”Spring constant”）、減衰値（”Damping constant”）、非伸長時のバネ長（”Unstretched spring length”）の３つです。
バネ定数としては通常は大きな値（剛なバネ）を与えます。接触状態になるとこのバネが有効になります。このバネ定数が小さすぎると接触してからも動く（すなわち、めり込む）ことになってしまいますが、大きすぎると数値計算の状態が悪くなります。このように接触をバネで表現しているため、接触後にどうしても振動してしまいます。減衰はそれを抑えるために与える必要があります。
非伸長時のバネ長は、その名の通り、バネがこれ以上縮まないという長さです。flange_a は flange_b にこれ以上近づけません。つまり接触状態になります。例えば、以下のモデルのように片側が固定壁であり非伸長時のバネ長が 0.3m であるとすれば、質量の右端と固定壁の距離が 0.3m 以下になれば接触状態ということになります。
ElastoGap のこの接触の仕組みは、ElastoGap コンポーネントの Modelica コードに以下のように組み込まれています。
contact = s_rel < s_rel0;
f_c = smooth(1, noEvent(if contact then -c*abs(s_rel - s_rel0)^n else 0));
f_d2 = if contact then d*v_rel else 0;
f_d = smooth(0, noEvent(if contact then (if f_d2 < f_c then f_c
else if f_d2 > -f_c then -f_c else f_d2) else 0));
f = f_c + f_d;
1 行目の s_rel は ElastoGap の伸縮量で両端の変位差（= flange_b.s – flange_a.s）であり、s_rel0 は入力である非伸長時のバネ長です。この s_rel が s_rel0 より小さい時、論理変数である contact が true （すなわち、接触状態）になります。
第 2 行では不連続状態で生じるイベントの処理のための関数である smooth や noEvent が使われていますが、理解しやすいようにそれらを除けば以下の式になります。
f_c = if contact then -c*abs(s_rel - s_rel0)^n else 0;
f_c はバネによる接触反力で、contact が true の時（接触状態の時）はバネの伸縮で生じる復元力（-c*abs(s_rel - s_rel0)^n）が代入されますが、非接触状態の時は 0 になります。ここで c は入力のバネ定数であり、n は補正係数です。ヘルプによると、金属球間の接触では n=1.5、金属板間では n=1.0 が推奨されています。
第 3、4 行では接触時の減衰力（f_d）を定義しています。これも、smooth と noEventを除くと以下のようになります。
f_d = if contact then (if f_d2 < f_c then f_c else if f_d2 > -f_c then -f_c else f_d2) else 0;
ここでは接触発生時の過大な減衰力の急激な発生を抑えるために、減衰力が復元力を超えないように定義しています。これにより、接触した瞬間に “引き剥がすような” 不自然な力が発生することを防いでいます。
MassWithStopAndFriction
MassWithStopAndFriction コンポーネントは Mass と同様の質量要素ですが、質量が動くことのできる範囲を限定することができ、かつ、質量の動きに対して静摩擦と動摩擦を定義することができます。ここでは可動域の制限に着目して紹介します（摩擦に関しては、次回で紹介する予定です）。
MassWithStopAndFriction では質量とその両側の “壁” を指定します。その質量は壁の内側では自由に動くことができます（ただし、摩擦を働かせることはできます）。また、質量は長さを持つことも可能です。これを指定するための MassWithStopAndFriction のパラメータダイアログは以下のようになっています。
まず、質量の右側と右側の壁の距離（”Right stop for (right end of) sliding mass” : smax）および質量の左側と左側の壁の距離（”Left stop for (left end of) sliding mass” : smin）を設定します。次は質量の長さ（” Length of component, from left flange to right flange” : L）と質量値（” Mass” : m）を指定します。長さはゼロでも構いません。例えば、Smax=7、smin=-7、L=6 だとすると。この質量は左右に ±4 m 移動できることになります。
“Velocity dependent friction” から “Exponential decay” までの入力は摩擦関連の係数です。摩擦を考慮しない場合は上図のように陽にゼロを指定しておきます。ブランクにしておくとディフォルト値が使われて摩擦が考慮されることになってしまいます。
MassWithStopAndFriction のこの接触の仕組みは、MassWithStopAndFriction コンポーネントの Modelica コードに以下のように組み込まれています。
stopped = if s <= smin + L/2 then -1 else if s >= smax - L/2 then +1 else 0;
when stopped <> 0 then
reinit(s, if stopped < 0 then smin + L/2 else smax - L/2);
reinit(v, 0);
end when;
1 行目では論理変数である stopped に接触状態を与えています。すなわち、左壁に接触した場合は -1、右壁への接触は +1、それ以外（自由に動いている間）はゼロです。
3 行目以降では stopped がゼロではなくなった瞬間に質量の変位（s）と速度（v）を初期化しています。すなわち、左壁への接触の場合は s は smin + L/2 に、右壁への接触の場合は smax - L/2 にそれぞれ設定されます。また、v は接触前の値に関わらず、ゼロに置き換えられます。
この when 文は if 文とは異なり、状態が変わって再度条件（stopped がゼロではなくなる）が満たされるまでは実行されません。
例１：ElastoGap と MassWithStopAndFriction の比較
紹介した二つのコンポーネントの接触の仕組みと相違を理解するために同じ結果になるような簡単なモデルを “無理やり” 作成してみました（添付：ElastoGap_MassStop_Comparison.scm）。
ここでは質量の長さはゼロとし、MassWithStopAndFriction の両側の壁までの距離（smin および smax）は共に 7.0m としているので、初期位置から両側に 7.0m 移動することが可能であることになります。
ElastoGap 側のモデルをこれに合わせるために、ElastoGap の Unstretched spring length（s_rel0）を -7 、s_rel の初期値を 0 とそれぞれ指定しています。これにより、ElastoGap の左右端は初期状態で同じ位置にあり（つまり、ElastoGap の初期長はゼロ）、flange_a は flange_b を超えて 7.0m 移動できることになります。ちなみに、左側の ElastoGap は移動量の正負が逆になるため、flange_a と flange_b の方向は右側の ElastoGap と同じになっています（鏡面対象にはなっていません）。
このモデルの、入力荷重としてサイン波を与えた場合の結果は以下のようになり二つのコンポーネントの挙動がほぼ一致していることが分かります。また、変位が両側とも 7.0m を “ほとんど“ 超えていないことが分かります。
ElastoGap の接触直後の “めり込み” や “振動” は ElastoGap で指定されたバネと減衰の値に依存しており調整の必要があります。例えば、ElastoGap の c と d をそれぞれ 10倍（c = 10000, d = 1000）にすれば以下のようにかなり抑えられることが分かります。
ちなみに、ここまで ”右” とか “左” とかを多用しており、また、見た目に理解しやすいように SlastoGap をモデル上は左右に配置していますが、これらは計算上は意味がありません。例えば、以下のようなモデルでも質量長さがゼロの場合は全く同じ挙動になることをご理解ください。
例２：ボールバウンシング（ボールの跳ね返り）
ElastoGap の例題として Demo ライブラリに含まれている BouncingBall モデルを紹介します。これは、ボールを空中で離し、重力によって落下して床に当たって跳ね返る、ということをシミュレーションしたモデルです。
モデルは以下の通りです。
Mass には重力を模擬した一定の力が与えられています。ElastoGap のダイアログは以下のようになっており、s_rel の初期値を 4.0、Unstretched spring length を 0.05 としていることから、4.0 の位置からボールを離し、壁から 0.05 離れた位置で跳ね返るように設定されていることが分かります（なぜ、0.0 ではなく 0.05 なのかは分かりません）。
結果は以下のようになります。
跳ね返る高さが徐々に低くなってくる効果は ElastoGap の減衰項（damping Constant）で表現されています。つまり減衰項をゼロにすればいつまでも跳ね続けます。
ところで、この例題の ElastoGap を MassWithStopAndFriction に置き換えることは可能でしょうか。結論から言うと、そのままでは使えません。仮想のバネとは言え Elastogap ではバネが働いているためボールが跳ね返されます。前述のように、MassWithStopAndFriction では、床に当たった瞬間にボールの速度がゼロに置き換えられてしまいます。跳ね返るには反発する前より少し（反発係数分）小さい速度を接触前とは逆向きに持つ必要があります。
このことから、MassWithStopAndFriction の Modelica コード内の reinit 文を変更すればよいことが分かります。すなわち、接触時には速度をゼロに初期化することになっている現在の reinit 文を以下のように変更します。
reinit(v, 0);
↓
reinit(v, -e*pre(v));
ここで e は反発係数で定数として値を入力します。また、pre() は直前のステップの値を返す関数です。これにより接触時には接触直前の速度とは反対向きの反発係数分小さくなった速度の値によって初期化されることになります。
このカスタマイズは MoCustomComponent ブロックを用いて行うことが可能で、以下のようなモデル（添付：BouncingBall_Cpmparison.scm）になります。
右側は ElastoGap による Bouncing Ball のモデルであり、左側がカスタマイズされた MoCustomComponent で置き換えたモデルです。反発係数を調整すると（e = 0.9264）以下のように左右の結果が一致します。
この Bouncing Ball のシミュレーションはいろいろな方法で試みられています。例えば、以下のサイトでは Modelica のライブラリを使わずに Modelica コードを直接書いて Boucing Ball の挙動を表現しています。
https://mbe.modelica.university/behavior/discrete/bouncing/
また、シグナルベースでモデル化した Boucing Ball が Activate のデモモデルとして以下の場所に含まれています。
"Activate" ⇒ "Activate Introduction" ⇒ "Boucing Balls" ⇒ "Bball"
Activate のアクティベーション機能を利用することにより、Boucing Ball のような不連続な挙動をシグナルベースモデルでも簡単に処理することができることがこの例題からお分りいただけると思います。
次回は、摩擦系のコンポーネントである、SupportFriction、Brake、MassWithStopAndFriction を紹介します。