物理シミュレーション目次

物理シミュレーション目次

 

AD

全方位に投げつけろ!

 昔、うんこ投げるゴリラいたな。

 いや、そうじゃない。

 

 前回の実装は運動方程式に基づいて動かしてるので、当然、上に投げると重力で落ちてくるようになってます。

 

 ただ、カメラポジションとか立方体の大きさなんかも、初期高度を示す変数altitudeに依存させちゃってるので、altitudeを0mとかにすると画面が変になります。なので、元々のaltitudeは画面最大高度ということにし、新しく高度用変数positionを追加することにしました。で、このpositionを逐次更新することにしdisplacementは廃止します。

	var altitude = 828;		//	初期高度画面最大高度 m
	var startDate;			//	開始時刻
	var lastDate;			//	最後の移動処理の時刻
	var velocity = 0;		//	初速度 m/s
	var displacement = 0;		//	変位 m
	var position = 0;		//	高度 m

 それと、空気抵抗も落ちるの前提だったので、これを進行方向に逆らう形に直す。

	function sign(value) {
		return (value < 0) ? -1 : 1;
	}

	function updateTime(deltasec) {		
		・・・
		var spv2 = s / 2 * p * velocity * velocity * sign(velocity);
		・・・
		velocity += (-9.8 + (-a * Lnv - b * spv2) / m) * deltasec;
		・・・
	}

 こうしないと、投げ上げる時に空気抵抗で上方向に加速しちゃうことになるんでな。

 新しく追加したsignメソッドは、引数で受け取ったvalueの値が正なら1、負なら-1を返すというもの。その中でやってる

		return (value < 0) ? -1 : 1;

 は三項演算子( ? : )を使ったもので、条件によって値を変えたい時に使うもんです。

 

    条件 ? 条件が真の時の値 : 条件が偽の時の値

 

 結果、以下のように書いたことと同じ意味になります。

	if (value < 0) {
		return -1;
	} else {
		return 1;
	}

 後は、displacementの削除、positionの追加に合わせた細かな調整なので、説明は省略。重要なのはupdateTime関数の計算の変更で、こうすることで速度が上向きなら下向きに、下向きなら上向きに空気抵抗が働くことになる。

 これで準備OKなんで、例えば20m/sの初速度を与えれば…

 

    var velocity = 20;                    //    初速度 m/s

 ↓こうなる

 

http://www.tetera.jp/xcc/ameba/free-throw.html

 

 なかなかいい動きします。

 で、ここまできたら、せっかくの3次元表示なんで、縦横斜め、八方になげられるようにしたいですよね。私はしたいです。

 

 全然関係無いけど、よく言われる四方八方の8や、八極拳の8は3次元空間のx,y,z3軸、左右、上下、前後の組み合わせの8だから、平面上で8方向想像してた人、間違いだからそれ。一極集中(ビッグクランチ)の太極拳と、全方位膨張(ビッグバン)の八極拳です。小学校の頃、二極拳〜七極拳とかもあるのかと、そんな風に考えていた時期が俺にもありました。

 
 なので高度を扱うだけの上下方向のみだった速度や加速度、位置を3次元に拡張します。

 スカラ:scalarからベクタ:vectorへ。

 スカラってのは、これまで使ってきた単独の量のことね。

 

  5とか、2とか

 

 ベクタってのは、複数のスカラで構成される量。

 例えば、x軸、y軸、2つの座標値をペアにした2次元ベクトルってのを習ったと思いますが、あれがベクタです。

 1つの点の位置を表現するのに2つのスカラ値(x、y)を使っていたでしょ。

 

  (5, 2)とか

 

 

 ベクトルもベクタも、英語のvectorをどう読むかであって、全く同じだけど、日本では、ベクトルは「向き」ってイメージが強いっすね。

 

 でも複素数をベクタで表現(実数部, 虚数部)したりもするんで、一概に「向き」って言っちゃうのはどうかと思われ。スカラ、ベクタ混在で構成するベクタとかもあるし。そもそもスカラ値のプラスマイナスからして数直線上の向きって言えなくもないし…

 ちなみに複素数の複素平面(ガウス平面というのじゃよ)を見て、あれ、これ2次元表現するのにちょうどよくね?って極座標使ってeのなんちゃらってやってsin、cosの式を簡潔に書いたりするとか、3次元でも使える複素数あるんじゃねとか言ってハミルトンが4元数とか言い出す話はまた今度。

 

 まあどっちでもいい。

 

 で、Three.jsにも3次元ベクタがオブジェクトとして定義されてるんで、こいつを使って、先のプログラムを書き直しました。

	var velocity = new THREE.Vector3( 0.7, 20, 0);	//	初速度 m/s
	var position = new THREE.Vector3( 0, 1.7, 0 );	//	初期位置 m

 THREE.Vector3はx,y,zプロパティを持つオブジェクトです。

 プロパティ?オブジェクト?って人は「走れJavaScriptちゃん!」から出直してこい!

 ま、それはいいとして、上のように書くと、velocityは(x、y、z)要素の順に(0.7, 10, 0)が設定されたTHREE.Vector3オブジェクトを示すようになり、positionは( 0, 1.7, 0 )が設定されたTHREE.Vector3オブジェクトを示すようになります。

 で、このpositionのx、y、z値をsphereのpositionにcopyメソッドを使って設定しています。実はTHREE.Meshの持つpositionプロパティもTHREE.Vector3なのでこんなことができる。「超高校級のMikuMikuDance」で話した「positionはカメラオブジェクトのプロパティで、3次元空間でのカメラ座標を示す。この3D座標を表現するプロパティ自身もオブジェクト」もTHREE.Vector3。

	function arrengeObjects() {
		・・・
		sphere.position.copy(position);

 それと床に立方体を使うのはやめて平面(THREE.PlaneGeometry)にしました。ここら辺は自分で調べてみてください。

 updateTimeメソッドもTHREE.Vector3に合わせて変更してるけど、x、y、zそれぞれの要素に対して運動方程式で計算してるだけで違いはないです。

 加速度のy成分だけ重力加速度が追加されるわけやね。

		var F = new THREE.Vector3();
		F.x = - nR * velocity.x 
			- D * (velocity.x * velocity.x) * sign(velocity.x);
		F.y = - nR * velocity.y 
			- D * (velocity.y * velocity.y) * sign(velocity.y); 
		F.z = - nR * velocity.z 
			- D * (velocity.z * velocity.z) * sign(velocity.z); 

		var acc = F.divideScalar(m);
		acc.y += -9.8;

 それと人間を投げあげるのは豪快すぎるので、バレーボール投げあげることにします。

 

  半径:10cm

  質量:200g

 

としました。それに合わせ、粘性抵抗

 

 

については、バレーボールにしたことだし素直にナビエストークの式使って

 

    F = 6πrηv

 

とします。結局人間用のαがわからんかったからね。

 

 THREE.Vector3のaddメソッドは、引数で渡したTHREE.Vector3のx,y,z値を加算するというものです。

 例えば、a、b、二つのTHREE.Vector3オブジェクトがあれば

 

        a.add(b)

 

でaのx,y,zは、それぞれbのx,y,zが加算された状態になります。

 

   a.x += b.x;

   a.y += b.y;

   a.z += b.z;

 

 こんな感じ。「+=」は左辺の変数(またはプロパティ)に右辺の値を加算という意味ね。

 multiplyScalarの場合は、引数で渡した値がx,y,z値にかけられます。

 

        a.multiplyScalar(2)

 

なら

 

        a.x *= 2;

        a.y *= 2;

        a.z *= 2;

 

とやったことと同じです。「*=」は「+=」の掛け算版です。

 でもってdivideScalarはその割り算版。

 cloneメソッドは複製で

 

    b = a.clone();

 

でbは、aのx,y,z値を持つ、新しいTHREE.Vector3オブジェクトを示すことになります。

 

 詳細はupdateTimeメソッド見てもらうとして、これでようやく斜め発射もできるようになる。

 で、ここまできたら床で跳ね返らせたいわけですよ。

 なので、床に当たると速度ベクタを運動エネルギー保存則使って変換します。

 

   v' = v - 2(v・n)n

 

   v':衝突後の速度ベクタ

   v:衝突時の速度ベクタ

   n:床の法線ベクタ

   

 v・nってのは内積です。

 THREE.Vector3にはdotメソッドとして用意されてます。便利。

 完全弾性衝突なんで、エネルギー減らず、減るのは空気中の移動による抵抗だけってことで、ボールはしつこく跳ね続けます。

 

 ↓それが、これだ!

http://www.tetera.jp/xcc/ameba/free-throw-3d.html

 

    ↓jsdo.it版

http://jsdo.it/reborn_xcc/ULWR

 

 完全弾性衝突とか内積とか、法線ベクタとか、詳細は次回!

AD
 
 

 JASRACが怖くて日和りました。

 Swiftの話だと思った?

 残念、ジュリーでした。

 

 前回(半年前だけどな)予告してた空気抵抗を考慮した自由落下を数値的に解くをやってみます。

 まずは、数字見ても面白くもなんともないので、cube-600.htmlを加工して、平面の上に球を落とすことにします。

 cube-600.htmlは「WebGLでいきますぜい」で紹介したページね。

 

 ↓これね

http://www.tetera.jp/xcc/ameba/cube-600.html

 

 何やってるかや、使い方は以下を順に読みましょう。

 

WebGLでいきますぜい

JavaScriptちゃんでプログラムしますぜい

 
 で、今回はこのcube-600.htmlを加工して、立方体(cube)に向けて、球を(sphere)を自由落下させてみます。
 だいたいドバイタワーくらいの高さ(828m)から落下させようと思うので、カメラポジションとかも調整が必要。高度を変えるたびに調整するのは面倒なので
 

 altitude(高度)

 

という変数を用意して、このaltitudeに設定された高さから、高度0mに配置したcubeまでを1画面にとらえるようにします。

 

 

 具体的には、カメラポジションとカメラの向きをaltitudeをもとに計算することで上のような画面になる。

    var altitude = 828;                    //    初期高度 m
   ・・・
    function threeStart() {
   ・・・
        camera.position.set(0, altitude * 1.1, altitude * 1.1);
        camera.up.set(0, 1, 0);
        camera.lookAt( {x:0, y:altitude * 0.55, z:0 } );

 ここで指定してる3次元座標の1はそのまま1mと考えてください。座標の単位をなんにするかは作る人の勝手です。m単位でもkm、mmなんでもいいけど、自分がどういった単位で座標を扱ってるかは、忘れないように。

 物理ではMKS単位系で計算式を表現することが多いのでmにしました。

 前回の抵抗を考慮した自由落下の式も単位はMKS単位系です。

 

 ↓こいつね

 

 

 cubeの大きさやsphereの大きさも微妙に調整しますた。

 ここら辺は、興味がある人は各自で

    function arrengeObjects()

を見てください。特に重要な話じゃないです。

 重要な点は、ページの表示から経過した時間によって、sphereを動かさなければいけないこと。

 なので、時刻を記憶するための変数を用意。

    var lastDate;                            //    最後の移動処理の時刻

 こいつに、ページ表示時に呼ばれるthreeStart関数内で

        lastDate = new Date();

とすることで、この処理が実行された時点の時刻を記憶させています。

 そして繰り返し呼び出されるloop関数で

        var now = new Date();
        var deltasec = (now.getTime() - lastDate.getTime()) /1000;

として、記憶した時刻(lastDate)からの経過時間を秒単位で計算します。getTime()はミリ秒単位なので1000で割ってる。

 これで経過時間が測定できたので、あとはこの経過時間を使ってオイラー法(前回説明した、移動距離を数値的に解くやり方の呼び名です)で位置を更新するわけですよ。そして最後に、次回のloop関数呼び出し時に、今の時刻からの差分を計算させるためにlastDateを更新する。

        lastDate = now;

 loop関数は、window.requestAnimationFrame(loop)によって、画面の更新間隔(大抵、毎秒60回)で呼び出されるんで、こんな面倒なことせずにdeltasecを1/60=0.17に固定してもいいんですけどね。実際、私のマシンではdeltasecは0.017秒くらいになってたし…

 ただ、毎秒60回の呼び出しはあくまで努力目標なのと、必ずしも画面の更新間隔(リフレッシュレートと呼ぶのだよ)が60Hz(ヘルツと呼ぶのだよ、いや、わかっとるわ)とも限らないんで、今回は真面目に上記のような経過時間測定をしてます。

 

 で、いよいよオイラー法。

 オイラー法で位置を更新する部分はupdateTime関数というのを作ってまとめています。

 当然、変数velocityとdisplacementは値を記憶するためにupdateTime関数外で定義ね。

    var velocity = 0;                 //    初速度 m/s
    var displacement = 0;             //    変位 m
   ・・・
    function updateTime(deltasec) {
   ・・・
        displacement += velocity * deltasec;
        sphere.position.set(0, altitude + displacement, 0);
    }

 velocityとdisplacementがオイラー法で計算される速度と変位です。

 altitudeの初期高度に現時刻での変位を加算して、現時刻で高度としてます。

 それをsphereのy座標に設定して位置を移動させている。

 

 で、ここで残念なお知らせ。

 粘性抵抗側を求める式

 
 
のαの適切な値が見つかりませんでした。
 どこ探しても球体での式

 

 6πrηv   π=3.14  r=球体の半径(m)

 

しか見つからない。Lの人間の体積はだいたい

 

 L = 0.06

 

みたいで、η(エータと読むのだよ)は1.8 x 10-5てのがわかったんだけどαが見つからない。

 

 ↓ここね

http://www.mterm-pro.com/machine-yougo/fluid-dynamics/water-air-bussei.html

 

 しょうがないので、球体の体積が0.06になるrを計算してα=76と決めました。

 

 

 もういっそのこと、粘性抵抗側は無かったことにしたいな。

 まあ、気を取り直して慣性抵抗側。

 

 

 こっちはスムーズに見つかって、βを形状を円柱とした1.2、ρ(ピーじゃないよローだよ)を1.2kg/m3と決定。

 

 ↓ここから

https://ja.wikipedia.org/wiki/抗力

 

 進行方向断面積は適当に0.33m2、質量mは60kgとして、いざオイラー法!

    function updateTime(deltasec) {        
        var a = 76;                          //        粘性側係数
        var b = 1.2;                         //        抵抗側係数(円柱)
        var L = 0.06;                        //        体積(人間、適当)    m3
        var n = 1.8 / 100000;                //        粘度(空気)  Pas
        var Lnv = L * n * velocity;
        var s = 0.33;                        //        進行方向断面積(人間、適当)  m2
        var p = 1.2;                         //        密度(空気)    kg/m3
        var spv2 = s / 2 * p * velocity * velocity;
        var m = 60;                          //        質量    kg
        velocity += - ((9.8 + (-a * Lnv - b * spv2) / m)) * deltasec;
        displacement += velocity * deltasec;    
        sphere.position.set(0, altitude + displacement, 0);
    }
 ドバイタワー(828m)から飛び降りたらどうなるかシミュレーションしてみた。
 
 ↓とうっ!
 ↓jsdo.itの方にも置いておく
 180km/hくらいで落下速度が安定しますな。
 てなわけで、ドバイタワーから飛び降りると19秒くらい落下を楽しめますよっと。
 空気抵抗ないなら、どんどん加速して460km/hくらいになって13秒で地上に激突。
 確かめたい人はupdateTime関数の変数a,bを0にすればいいと思うよ。
 実際、成層圏から飛び降りた映像でも、空気の薄いところでは音速超えて、そこから落ちていって220km/hくらいで落ち着いてるのがわかる。
 以外に遅い。
 

 

 

 「途中でパラシュート開いたらどうなるの?」てのを試すなら、高度が100m切ったくらいのところでsの値を20くらいにするといいでしょう。
 一気に20km/hまで減速して、無事地上に降りられます。
 ま、パラシュートが一瞬で開くわけでもないから、0.5秒くらいで全開になるようにしても面白いかと。そこらへんは各自で。
 

 画面下にある現在の速度や経過時間、高度を更新する

    function updateInfo(now)

に関しても説明はパスします。計算結果を画面に表示してるだけです。

 「走れJavaScriptちゃん!」を読むとか、JavaScript getElementById innerHTMLあたりで検索すればすぐわかります。

 

 あとthreeStart関数でやってる
        light = new THREE.DirectionalLight(0xFFFFFF, 100);
        light.position.set( 0, 2000, 0 );
        scene.add(light);
   ・・・
        renderer.setClearColor(0xEEEEEE, 1.0);
        renderer.shadowMap.enabled = true;
        light.castShadow = true;
        light.shadow.camera.left = -altitude;
        light.shadow.camera.bottom = -altitude;
        light.shadow.camera.right = altitude;
        light.shadow.camera.top = altitude;
        light.shadow.camera.far = 3000;
        cube.receiveShadow = true;
        sphere.castShadow = true;
という部分も、画面上で影をつけるための処理で単なる見栄え。
 気になる人は直接Threeのサイトで調べるか、Three.js 影 あたりでググりましょう。
 
 ↓Threeのサイト
AD