２つの正電荷をおびた粒子がバネでつながっている物体のシュミレーション

こんにちは。
２つの正電荷をおびた粒子がバネでつながっている物体について、粒子をランダムな位置において離し、停止（ダンピングあり）するまでの過程をシミュレーションしたいのですが、どうもうまく行きません。

ちなみに空間は二次元（平面）です。

ある資料をもとに下記のVBAをコーディングしたのですが（簡単にしてあります）、力の向きをどういう風に反映させたらいいのかわかりません。

どうしたらいいでしょうか。。。　漠然としていてすみません。

粒子＝Node1, Node2

{

'クーロン反発力を計算するための関数
Public Function AddCoulombRepulsion()
Dim xyIdx As Integer
Dim distance As Double

For xyIdx = 1 To 2
'距離を計測
distance = Abs(Node1.position(xyIdx) - Node2.position(xyIdx))
'クーロン反発力を計算
Node1.netForce(xyIdx) = Coulomb定数 * Node1.電荷 * Node2.電荷 / distance ^ 2
Next xyIdx

End Function


'フック張力を計算するための関数
Public Function AddHookeAttraction()
Dim xyIdx As Integer
Dim distance As Double

For xyIdx = 1 To 2
'距離を計測
distance = Abs(Node1.position(xyIdx) - Node2.position(xyIdx))
'バネによる張力を計算
Node1.netForce(xyIdx) = ( -1 ) * (Spring1.バネ定数 * distance)
Next xyIdx

End Function

}
上記をNode1とNode2を逆にして再度行います。

{
'netForceから速度を計算
For xyIdx = 1 To 2
Node1.velocity(xyIdx) = (Node1.velocity(xyIdx) + timestep * Node1.netForce(xyIdx) / Node1.Mass) * Damping '※0 < Damping < 1
Next xyIdx

'velocityから位置を計算
For xyIdx = 1 To 2
Node1.position(xyIdx) = Node1.position(xyIdx) + timestep * Node1.velocity(xyIdx)
Next xyIdx

'totalKineticEnergyに加算
For xyIdx = 1 To 2
totalKineticEnergy = totalKineticEnergy + Node1.Mass * (Node1.velocity(xyIdx)) ^ 2

Next xyIdx
}
Node1をNode2にして上記を繰り返します。


これを、totalKineticEnergyが一定値より小さくなるまでループさせる予定です。(totalKineticEnergyは毎ループの頭で0に初期化します）停止したときのNode1の位置とNode2の位置が知りたいのですが。。。


ご助言よろしくお願いいたしますm(_ _)m

No.3 ベストアンサー 回答者： htms42 回答日時： 2008/09/09 09:55

文章を読む限り、バネは1つですね。

バネが振動しながら減衰していく途中の時間変化を数値的に求めようとしているのだと思います。でもこれはシミュレーションとは言わないと思います。ランダムというのは任意の初期条件という意味ですね。確率的に起こる現象を想定しているわけではないでしょう。
バネでぶら下がった錘が振動しているのとの違いは重力かクーロン力かです。重力は距離によらず一定ですがクーロン力は距離と共に小さくなります。つりあいの位置に関して力が非対称ですね。
振動の特徴を知ることと減衰項の考慮の仕方を勉強するためであれば
１次元で（ぶら下がっているバネ、水平に置かれているバネ）でやられてもいいのではと思います。
２次元ということは特に意味を持っていないという印象です。

速度に比例する抵抗がある場合の自由落下は微分方程式を解くことができますので数値的に解く方法の確認に使うことが出来ます。
timestepごとに力、加速度、速度、位置と求めていけばいいだけです。
この場合は振動ではありませんからダンピング項の扱いの練習になります。
運動方程式は　ｍａ＝ｍｇ－γｖ
ｔ＝０でｖ＝０　として解くと　ｖ＝(ｍｇ/γ)(１－exp(－γｔ/ｍ))　です。

減衰項のある場合の調和振動も微分方程式が解けているはずですから対応が可能でしょう。この場合は振動ですから停止の判断が問題になります。「運動エネルギー＝０」という条件では不足です。振動の途中で速度＝０が起こるからです。
「全エネルギー＝弾性エネルギー＋運動エネルギー」
です。止まった時はつりあいの位置になっていますから
「全エネルギー≧つりあいの位置での弾性エネルギー」
が静止の判定基準になるはずです。
＞totalKineticEnergy = totalKineticEnergy + Node1.Mass * (Node1.velocity(xyIdx)) ^ 2
この式は運動エネルギーだけしか考慮されていませんから停止の判定条件としては不足だと思います。

この回答への補足

本当は、扱っているのは複数の粒子が複数のバネでつながったネットワークです。
現実世界に即していなくても、ネットワークがうまいこと広がればそれでいいのですが。。。

補足日時：2008/09/10 22:41

No.2 回答者： rabbit_cat 回答日時： 2008/09/08 23:38

＞一応(xyIdx)が1または2で、それぞれxとyを表しています。

なるほど。じゃあ、このままでＯＫなのでは。

ただし、
＞Dampingは0次元？
＞だと思います（＾＾；）
＞Node1.velocity(xyIdx) = (Node1.velocity(xyIdx) + timestep * Node1.netForce(xyIdx) / Node1.Mass) * Damping
この式だとDampingは無次元になってないですよ。（パッと見は無次元になっているように見えますが、よく考えるとおかしい）
この式だと、ループ毎に、速度がDamping倍になっていくわけですが、ループ１回の時間がtimestep秒なわけですよね。ここで、もしtimestepが1/2になったら、Dampingは1/2も1/2にならないとおかしくないですか？

そもそも、ダンピングというのは、普通は、速度に比例する力が物体にかかる場合をいうことが多いのですが。
Node1.velocity(xyIdx) = Node1.velocity(xyIdx) + timestep * (Node1.netForce(xyIdx) - Damping*Node1.velocity(xyIdx)) / Node1.Mass
とか。
示された式だと、速度自体が、そのときの速度に比例して減っていくことになっています。本当にそういう系ですか？

この回答への補足

う～ん、どうなんでしょう(^^;)
元ネタにあった通りにやってみたのですが。。。

元ネタはこれです↓(wikipedia: http://en.wikipedia.org/wiki/Force-based_algorit …

Each node has x,y position and dx,dy velocity and mass m. There is usually a spring constant, s, and damping: 0 < damping < 1. The force toward and away from nodes is calculated according to Hooke's Law and Coulomb's law or similar as discussed above.

set up initial node velocities to (0,0)
set up initial node positions randomly // make sure no 2 nodes are in exactly the same position
loop
total_kinetic_energy := 0 // running sum of total kinetic energy over all particles
for each node
net-force := (0, 0) // running sum of total force on this particular node

for each other node
net-force := net-force + Coulomb_repulsion( this_node, other_node )
next node

for each spring connected to this node
net-force := net-force + Hooke_attraction( this_node, spring )
next spring

// without damping, it moves forever
this_node.velocity := (this_node.velocity + timestep * net-force) * damping
this_node.position := this_node.position + timestep * this_node.velocity
total_kinetic_energy := total_kinetic_energy + this_node.mass * (this_node.velocity)^2
next node
until total_kinetic_energy is less than some small number //the simulation has stopped moving

補足日時：2008/09/09 06:52

No.1 回答者： rabbit_cat 回答日時： 2008/09/08 00:00

＞Node1.velocity(xyIdx) = (Node1.velocity(xyIdx) + timestep * Node1.netForce(xyIdx) / Node1.Mass) * Damping

Damping という定数の（物理量としての）次元はなんですか？
timestepがかかってないとかなり変な感じです。

＞力の向きをどういう風に反映させたらいいのかわかりません
単純にやるなら、粒子の位置や力を２次元ベクトルで表わせばいいです。

この回答への補足

Dampingは0次元？
だと思います（＾＾；）

一応(xyIdx)が1または2で、それぞれxとyを表しています。

補足日時：2008/09/08 05:51