13.滑り・転がり

拘束された剛体の運動方程式は、式(8)で与えられる。運動方程式を確定させるには拘束条件が必要であるが、曲面上を滑る場合は式(18)であり、転がる場合は式(19)となる。
この章のシミュレーション(クリックで計算開始)

運動方程式(1)の拘束力 fc と拘束トルク τcg を求めたい

床の上に置かれた剛体の運動を計算したい。剛体には、摩擦なく滑る(=滑り拘束)、あるいは、全く滑らずに転がる(=転がり拘束)という拘束条件が課されているとする。前章では、剛体上の点を固定することを考えたが、その場合とはかなり違った運動になりそうである。
運動方程式を求めるための戦略としては、2種類が考えられる。1つ目は、前章のように自由な速度を求め、それに対する運動方程式を立てる方法。2つ目は、自由な剛体の運動方程式(第11章の11.1節)に、拘束力の総和 fc と拘束トルク τcg を加えたもの:(添え字 g は剛体の重心を基準にしていることを表す) m¨xg=f+fcddt(Igω)=τg+τcg(1) を求める方法である。 fc,τcg を求めれば、運動方程式が確定することになる。
今考えている状況では、床と接触している点が時々刻々変わるので、1つ目の方法は難しそうである。よって、2つ目の、式(1)を用いる方法を採用しよう。
この章では、まず、剛体に拘束を追加した時の運動方程式を一般的に求める。その後、床と剛体が接しているための条件を考えることにより、滑り・転がり拘束の条件式を具体的に書下す。これらを、以下のように4つの節に分けて議論する。 13.1:13.2:13.3:13.4:

13.1 拘束が追加された剛体の運動方程式

計算を見やすくするため、式(1)を以下のようにまとめておく: ddt(μΩ)=Φ+Φc(2) ただし、各々の記号は以下のように定義している: μ[m00Ig],Ω[˙xgω],Φ[fτg],Φc[fcτcg] 運動方程式を決めるためには、 Φc (これも拘束力と呼ぶことにする)を求めればよい。(重心を基準にしているので μ などは μg のように書いたほうが誤解はないが、見づらいので添え字 g は省略する。)
この節では、剛体に拘束条件を課した時の拘束力 Φc がどうなるかを一般論を考え、運動方程式が式(8)で与えられることを見る。そうすれば後の節で、滑り・転がりの場合の拘束条件を求めれば、それらの運動方程式が確定することになる。

ダランベールの原理(4)により、運動方程式は式(8)

剛体の自由な速度 Ω に対して、以下のような拘束条件がかかっているとする:B は行列) BΩ+b=0(3) これはちょうど、 G(t,X)=0 という形の拘束条件の両辺を t で微分して、速度 ˙X に対する拘束条件にしたもの GX˙X+Gt=0 と同じ形である。
式(3)から、ダランベールの原理(第8章の8.2節)に対応する条件を求めれば、拘束力 Φc の向きが決まるはずである。まず、剛体の運動エネルギー E を、 Ω を用いて書くと、以下のようになる:(以下の【13.1-注1】) E=12ΩTμΩ 第8章での議論と同様に、拘束力 ΦcE を変化させないので ΩTΦc=0 となる。これが BΩ=0 (式(3)において b=0 としたもの)を満たす全ての Ω に対して成り立たなければならない: ΩTΦc=0BΩ=0 即ち、 Φc は、未定乗数 λ を用いて Φc=BTλ(4) と書ける。これが、剛体に拘束条件を追加したときのダランベールの原理である。
よって、運動方程式(2)に式(4)を代入すれば ddtμΩ=Φ+BTλ(5) となる。未知数は λ のみである。 λ を求める方法も第8章の場合と同じである。即ち、式(5)と「拘束条件(3)を時間微分したもの」を連立して、 ˙Ω を消去すればよい。実際に計算すると、求める運動方程式は、以下の【13.1-注2】の式(8)となる。

【13.1-注1】剛体の運動エネルギー

剛体の運動エネルギー E は、以下のようになる: E=12ΩTμΩ=12m˙xg2+12ωTIgω(6) (最後の式では、第1項が重心運動のエネルギー、第2項が回転運動のエネルギーという形にうまく分離している。)

導出

剛体を構成する質点要素の速度 ˙xiΩ の関係式は以下である:(第11章の11.1節) ⎢ ⎢˙x1˙x2⎥ ⎥=⎢ ⎢1x1g×1x2g×⎥ ⎥ΩAΩ これを使って、運動エネルギー E の中の ˙xiΩ で表せばよい: E=12imi|˙xi|2=12[˙xT1˙xT2]⎢ ⎢m1m2⎥ ⎥⎢ ⎢˙x1˙x2⎥ ⎥=12ΩTAT⎢ ⎢m1m2⎥ ⎥A=μΩ 最後の式の μ の定義は、11.1節で述べたものである。

【13.1-注2】拘束された剛体の運動方程式:式(8)

剛体に、拘束条件 BΩ+b=0(7) がかかっている時、運動方程式は以下のようになる: ddt(μΩ)=Φ+BTλλ(Bμ1BT)1[˙BΩ+˙b+Bμ1(Φ˙μΩ)](8)(9) なお、 ˙μ は以下のように書ける:(第11章の【11.1-注2】) ˙μ=[000ω×IgIgω×]
これを計算するには、初期値が式(7)を満たしている必要がある。しかし、式(7)が可積分の場合には、積分したものも満たさなければならない(以下の【13.1-注3】参照)

導出

式(8)の未定乗数 λ が、式(9)になることを示せばよい。そのためには、式(8)左辺の時間微分を実行したもの: ˙μΩ+μ˙Ω=Φ+BTλ˙Ω=μ1(Φ+BTλ˙μΩ) を、拘束条件(7)の時間微分: ˙BΩ+B˙Ω+˙b=0 に代入して、 λ=[] の形にすればよい。そうすれば式(9)に一致する。

初期値に対する拘束条件

運動方程式(8)が得られたわけだが、これを解くためには、拘束条件を満たすように初期値を取っておく必要がある。そのように初期値をとっておけば、運動方程式(8)の解は、任意の時刻において、自動的に拘束条件を満たし続ける。
初期速度 Ω に対する拘束条件は、式(7)である。しかし、拘束条件は Ω だけに課せられるわけではなく、剛体の位置 xg と姿勢 θ に対しても拘束条件が課せられているはずであり、初期位置・初期姿勢はこれを満たさなければならない(例えば床に接しているという条件)
もし、拘束条件が、速度を含まない位置 xg と姿勢 θ の関係式 G(t,xg,θ)=0 の形にかけたならば、初期値が満たすべき拘束条件は、第8章のように、 G=˙G=0 である。その場合、拘束条件(7)は ˙G=0 に対応している。問題は、対応する G が常に存在するとは限らないということである。存在しない場合、その拘束条件は非可積分であるという(以下の【13.1-注3】)
非可積分な拘束条件の場合、初期値に対する拘束条件は、式(7)のみとなる(あるいは部分的に非可積分という場合もある)。従って、式(7)の形で拘束条件が与えられている場合、正しい初期値を与えるためには、その式が可積分かどうかを判定しなければならない。(【13.1-注3】の補足で述べた様に、この章で扱う問題の範囲では、数学的な議論なしで判定できる。)

【13.1-注3】非可積分な拘束条件

一般に、拘束条件 B˙X+b=0(10) が、何らかの拘束条件 G(t,X)=0 の時間微分 ˙G=0 で表せる場合[1]、式(10)は可積分(あるいはホロノミック)であると言う。逆に、表せない場合、非可積分(あるいは非ホロノミック)であるという。

補足

[1] ˙G=0 と表せたとしても、式(10)の左辺と ˙G が一致するとは言えない。というのも、式(10)に両辺に適当な関数を掛けてもよいからである。従って、 ˙G=0 で表せないことを示すためには、どのような関数を掛けたとしても、式(10)の左辺が ˙G の形にかけないことを言う必要がある。(第10章の【10.3-注1】で、角速度 ω に対し、 ˙θ=ω となるような回転の自由な座標 θ が存在しないことを示した。その時は偏微分の可換性を課すだけで示せたが、今の場合、関数をかける自由度がある分、より複雑である。)
以前の章で出てきた拘束条件は、最初から G が分かっている可積分なもののみであった。しかし、この章で扱う転がり拘束は、非可積分である。例えば、平坦な床の上にボールを置いて転がすことを考える。この時、床との接触点 xc を保ったままボールの姿勢を変えることは、自由にはできない(鉛直方向を軸とした回転のみが可能)。これは、「滑らない」という拘束条件のためである。しかし、「ボールを転がして移動させてから元の位置 xc に戻す」という操作と合わせると、実は、ボールの姿勢を任意に調節することができる(実際にボールを床の上で転がしてみれば分かる)。これが意味することは、 xc を決めても、剛体の姿勢は決まらないということである。即ち、剛体の位置・姿勢に対する拘束条件は、床に接しているという条件のみであり、「滑らない」という追加の拘束条件は、速度に対する拘束条件であって、位置・姿勢に対しては拘束を与えないのである。

13.2 剛体と床が接しているための条件

この節では、剛体と床が接触し続けるために、剛体の速度 ˙xg,ω が満たすべき拘束条件が、式(14)で与えられることを見る。滑り拘束と転がり拘束は両方とも、この条件を満たさなければならない。簡単のため、床と剛体は1点のみで接しているとする。

モデリング:時刻 t での剛体の形状は式(11)

床の形状を G(t,x)0 で表すことにする。即ち、これを満たす x の集合が床を構成する。 t に依存しているように書いているのは、時間とともに床が変形してもよいことを表している。また、モデル位置における剛体の形状を ˜H(x)0 で表すことにする。即ち、モデル位置における全ての質点要素が、この不等式を満たす。モデル位置は任意であり、床と接触するようにとる必要はない。
時刻 t において、剛体の位置および向きが xg,R であるとするxg は重心位置、 R は回転行列)。この時の剛体の形状 H(t,x)0 は、モデル位置での形状 ˜H を用いて、以下のように書ける:˜xg はモデル位置での剛体の重心) H(t,x)˜H(RT(xxg)+˜xg)0(11) これを示すには、モデル位置での質点要素 ˜xi が満たす式 ˜H(˜xi)0 に、時刻 t での質点要素の位置 xi=xg+R(˜xi˜xg) を代入して、 ˜xi を消去すればよい。 H が時間に依存しているように書かれているのは、剛体の運動によって、位置や向きが変わることによるものであり、剛体自体が変形するわけではない。

xg,R に対する拘束条件は式(12)、 ˙xg,ω に対する拘束条件は式(13)

剛体と床が接触しているための条件を考える。剛体と床の接触点を xc とおく(接触は1点のみと仮定する)xc が満たす条件は、「 xc が剛体と床の両方の表面にあり、かつその点での接平面が一致する」ことなので、以下のように書ける:ˆGG の大きさを正規化したもの) G=0H=0ˆG+ˆH=0(12) G,H はそれぞれ G,H が大きくなる方向を向くので、右図のように互いに逆を向くことに注意。) G,H は、ともに (t,xc) での値である。上式を満たす xc が存在するように H が選ばれている必要があるので、この式は、式(11)を通して剛体の配置 xg,R に対する拘束条件にもなっている。
剛体と床が接触し続ける時、上式は任意の時刻の xc(t) に対して成立するので、同式の時間微分も成立しなければならない: Gt+(G)T˙xc=0Ht+(H)T˙xc=0(ˆGt+ˆHt)+(ˆGx+ˆHx)˙xc=0⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪(13) H の時間微分には ˙xg,ω が含まれているので、この式は、 ˙xc に対する条件であると同時に ˙xg,ω に対する拘束条件にもなっている(上述のように、 H の時間依存性は、 xg,R の時間依存性によるものである)

式(13)の分離→ ˙xg,ω に対する拘束条件:式(14)

ところで、剛体の位置・向きとその速度は xg,˙xg,R,ω を与えれば決まるので、これらを与えれば、 δt 秒後の剛体の位置も決まる。よって、接触点の移動量 δxc も決まるので、接触点の速度 ˙xc も決まる。よって、式(13)は、「 ˙xg,ω に対する拘束条件」だけでなく「 xg,˙xg,R,ω から ˙xc を決定するための方程式」も含んでいるはずである。両者を分離したい。
まず、式(13)の第1式と第2式を、式(12)の第3式と連立すると、 ˙xc を含む項が消せて以下のようになる: 1|G|Gt=1|H|Ht(14) ˙xc がなくなったので、これが ˙xg,ω に対する拘束条件になっている。

式(13)の分離→接触点の速度 ˙xc :式(15)

次に、接触点の速度 ˙xc についてである。拘束条件(13)の第3式は3成分の式なので、 ˙xc=[] の形に変形できるだろうか。実はできない。実際、同式の全ての項に単位ベクトルの微分が現れているが、これらは元の単位ベクトルと直交する(第8章の【8.3-注1】)。そのため、同式に (G)T (またはそれに平行な (H)Tを左乗すると自動的に 0 となる。よって、3成分全てが独立というわけではない。
実際に ˙xc を求めるには、拘束条件(13)の第3式と、式(13)の別の1つの式と合わせればよい。例えば、第1式に G をかけた G(G)T˙xc=GtG を同第3式に辺々加えると、 ˙xc にかかる行列は可逆になるので ˙xc が求まる: ˙xc=K1(ˆGt+ˆHt+GtG)KˆGx+ˆHx+G(G)T(15)
以上により、式(13)の4つの条件は、拘束条件(14)と接触点の速度を与える式(15)に、過不足なく分離することができた。

13.3 滑り拘束

この節では、滑り拘束に対する拘束条件が式(18)となることを見る。これにより、すでに求めた運動方程式(8)が確定する。それを用いて、具体的な計算を行う。

拘束条件(14)を式(7)の形にする→式(18)

滑り拘束の拘束条件は、剛体と床が接触しているという前節の条件(14):(再掲) 1|H|Ht1|G|Gt=0(16) そのものである。よって、運動方程式(8)を決めるには、これを式(7)の形BΩ+b=0に変形すればよい。床が静止している場合、直感的には、「接触点 xc にある剛体の質点要素の速度が、拘束面に平行」という形になりそうである(そうでなかったら、めり込んだり離れたりしてしまう)。
その変形を行うには、式(16)の Ht の中に含まれている ˙xg,ω を表に出す必要がある。そのためには、 H˜H の関係式(11):(再掲) H(t,x)˜H(RT(xxg)+˜xg) を使えばよい。この式の両辺の全微分を取ると、 ˜x=RT(xxg)+˜xg とおいて Htδt+Hxδx˜H˜xδ˜x˜H˜x(˜xtδt+˜xxRTδx) となる。両辺の δt,δx の係数がそれぞれ一致するので以下を得る: Ht=˜H˜x˜xtHx=˜H˜xRT 後は、第2式を使って、第1式の緑字部分(要するに ˜Hを消去し、計算を進めればよい:tR,xg だけに作用することに注意) Ht=HxR˜xt‖ ‖ ‖ ‖ ‖ ‖ ‖ ‖ ‖ ‖ ‖˜xt=t[RT(xxg)+˜xg]=˙RT(xxg)+RT(0˙xg)+0=RT[ω×(xxg)+˙xg]=(H)T[˙xg(xxg)×ω](17)
拘束条件(16)は x=xc の時に成立する式なので、式(17)において x=xc としたものを式(16)を代入すればよい。その後、 Ω をくくりだすと、拘束条件は以下のように書ける: BΩ+b=0B(G)T[1(xcxg)×]bGt⎪ ⎪ ⎪⎪ ⎪ ⎪(18) (式(12)の第3式 ˆG+ˆH=0 を使って HG で置き換えた。)この B,b を用いることにより、運動方程式(8)が確定する。なお、床が静止している場合、 Gt=0 なので、この式(18)は、接触点の速度が拘束面と平行になっていることを意味しており、上で推測した通りになっている。

計算方法

以上をまとめると、滑り拘束の場合には以下のように計算すればよい:
  • まず、床の形状 G(t,x) と、モデル位置での剛体の形状 ˜H(x) を定義する。
  • 次に、拘束条件(12)を満たす接触点 xc が存在するように、 t=0 での初期位置 xg,R をとる。初期速度 ˙xg,Ω については、拘束条件(18)を満たすようにをとる。
  • その後の運動は、運動方程式(8)から計算できる。ただし、 B,b およびその微分は、 xc における値である(微分すると ˙xc が出てくるが、式(15)から求められる)
なお、接触点 xc は、 xg,R から直接計算することも原理的には可能であるが、ステップ毎に直接計算するのは、一般に困難である。その場合、式(15)の ˙xc を用いて xc の時間変化を同時に計算していくこともできる: xc(t+δt)xc+˙xcδt
前述のように、今考えているのは、常に1点で接触している状況である。接触点が複数ある場合、あるいは運動の途中で増える場合(衝突運動になることが多いだろう)を扱うことはできない。逆に、接触点が減る場合も扱えない。例えば、本来であれば空中に飛び上がってしまうような場合でも、床から離れないようにする拘束力が働くため、床にくっついたままになる(拘束力の総和 fc を計算しておいて、 fc の向きが床にひきつける方向になった瞬間に、 Φc0 にして自由な運動に切り替えれば対応できる)。

計算例

例として、球を x,y,z 方向からつぶした楕円体の剛体 ˜H(x)x2a2+y2b2+z2c21 を取りあげるa,b,c は定数)。密度は一様とする。モデル位置での慣性モーメント ˜Ig は、剛体を質点要素に分解して数値的に近似計算してもよいが、今の場合には解析的に計算することができ、以下のようになる: ˜Ig=m5b2+c2c2+a2a2+b2 m は剛体の質量である。導出は第14章で行う。
数値計算を行うと下図のようになる。最も長い軸の半径が 0.3m となるようにしている。

13.4 転がり拘束

この節では、転がり拘束での拘束条件が、式(19)となることを見る。前章の滑り拘束に、滑らないという拘束条件を加えればよい。

拘束条件は式(19)

滑らないということは、接触点 xc において剛体と床が相対速度を持っていない、即ち、「 xc に位置している剛体上の質点要素 xc,H 」の速度 ˙xc,H=˙xg+ω×(xcxg) と、「 xc に位置している拘束面上のxc,G 」の速度 ˙xc,G が等しいということである: ˙xg+ω×(xcxg)=˙xc,G よって、式(6)の形に変形すると BΩ+b=0B[1(xcxg)×]b˙xc,G(19) となる。
なお、 ˙xc,G を求めるためには、床の上の各点の速度を与える必要がある。しかし、動く壁との衝突(4.2節)の時にも述べたが、床の表面を定義する G=0 という式の中に、床の水平方向の速度の情報を含めることはできない。そのため、床が動いている場合には、追加で速度の情報を与えてやる必要がある。一方、床が静止している場合には、単に ˙xc,G=0 とすればよい。
式(19)の B は、滑り拘束の条件(18)の B を含んでいる。よって、転がり拘束の場合の計算方法は、滑り拘束の場合の拘束条件(18)を、上式(19)で置き換えるだけである。
剛体の初期位置・初期姿勢に対する拘束条件は、床と剛体が接していることのみなので、滑り拘束の場合と同じである。従って、初期位置に対する拘束条件よりも、初期速度に対する拘束条件(19)のほうが数が多い。即ち、転がり拘束は非可積分(【13.1-注3】)である。

計算例

床が静止している場合、数値計算を行うと下図のようになる。剛体の設定は、滑り拘束の場合と同じである。