というか結論から言いますと、やべぇwシミュレーションぱねぇわwって感じでした。
見た目は限りなく地味ですが、実践できた運動と同一の運動を解析的に解こうとすると(つまり方程式を解いてバシッと運動を決めようとすると)
どれほど面倒なことかということは日々の勉強から知っています。
だからいとも簡単に(コード上はそこまでですが、少なくとも厄介な微分方程式を解くよりははるかに簡単に)運動が再現できたときは驚きました。 実行ファイルでは順に
単純なばね振り子(ただし、計算上は平面上の振動として行っています。おんなじことなので)
連成振動(質点数2ただしいくらでも増やせる)
よくわからん複雑な運動(おそらく解析的に解くのはしんどい)
となってます
さて、どうやったかですが、やったことは高校で習うばねの運動そのものです。
まずばねの性質は以下のようなものです
ばねの力(N) = ばね定数(N/m)×ばねの伸び(m)
つまり、ばねは伸びれば伸びるほど、それに比例して力が強くなるというものです。
ばねは以下のように宣言されています
class Spring:public Entity{
public:
Spring(
Physics::Phy k, Physics::Phy l,
std::shared_ptr e1, const Linear::Vector>& s1,
std::shared_ptr e2, const Linear::Vector>& s2
)
:Entity(
Linear::Vector>(),
Physics::Phy(0.0),
Linear::Vector >()
),
k_(k),l_(l),
e1_(e1), s1_(s1),
e2_(e2), s2_(s2),
circleNum_(11), r_(1)
{}
void AddForce(const Linear::Vector> &f){}
void Update();
void Draw();
private:
Physics::Phy k_;
Physics::Phy l_;
std::weak_ptr e1_;
const Linear::Vector>& s1_;
std::weak_ptr e2_;
const Linear::Vector>& s2_;
int circleNum_;
int r_;
};
ばね定数と、自然長(伸び0のときの長さ)
およびばねの接続先二つを与えます(それぞれ接続先のオブジェクトと、そのオブジェクトのどの位置を参照するのかを与えます)
AddForceは基底クラスEntityの仮想関数で、外力が与えられたときの振る舞いを記述しますが、ばねには特にすることがないので何も書きません。
以下は毎フレーム呼ばれるUpdate関数の定義です(UpdateはEntityの純粋仮想関数です)
void Spring::Update(){
shared_ptr p1, p2;
if(! ((p1 = e1_.lock()) && (p2 = e2_.lock())) ) return;
Vector> d = s2_- s1_;
Phy x = d.Length();
Phy dx = x - l_;
d.Normalize();
Phy f = dx*k_;
p2->AddForce(-d.GetValue()*f);
p1->AddForce(d.GetValue()*f);
}
上2行は接続先のオブジェクトがまだ生きてるかを調べているコードで、指して意味はありません(生きていればp1,p2にポインタが代入されます)
3行目からですが、まず二つの接続されている位置の相対位置ベクトルを求めます
そしてその長さ(つまりばねの端から端の距離)を求め
自然長の長さを引くことで伸びを得ます
また相対位置ベクトルは正規化(つまり長さ1)にして、これをあとで力の方向に変えられるようにします
そして、ばね定数×伸びで力を導出した後
二つの物体に力を加えるというわけです
d.GetValue()はdの無次元量を返します(d*fでは次元が変わるので計算できない)
これで、ばねの運動が再現できるようになりました。
ただし、この実装ではあくまでバネは質量0で、かつ接続されている点同士を直線で結ぶだけという簡単な条件です。
梁の曲げとかを行うにはもうすこし工夫が要りそうな気がします
次は拘束条件でも考えてみようかな?
今回のプロジェクト