研究メモ ~python2Dシミュレーション編~
姿勢拘束関数におけるバネダンパ係数をなんとなく決めるために、pythonでアメンボジャンプを作ってみた
（シミュレーションするためにシミュレーションしているという謎は一旦置いておく）
？？？？？
アクロバティックすぎる
たぶんこれはバネ係数がデカすぎ
左右で物性値も運動方程式も対象的なはずなのに動きが違うのはどうしてだ？
昔のうまくいったファイルと見比べてみたら、maxInputのあるなしで違いがあったので、一度コメントアウトしたらちゃんと対称に動いた
よく見てみると、
if (t_femurR > ctrl.maxInput): t_femurR = ctrl.maxInput
if (t_femurL > ctrl.maxInput): t_femurL = ctrl.maxInput
これではRもLも同じ回転方向になるじゃん
修正
if (t_femurR < -ctrl.maxInput): t_femurR = -ctrl.maxInput
if (t_femurL > ctrl.maxInput): t_femurL = ctrl.maxInput
うーん絶対値の方がよさそうだが、一旦これで
OpenFOAMの方ではロドリゲスの回転行列を使っているため、正方向が剛体によって別
一方こちらでは角度もトルクも向きが共通なので（どうせ2方向しかないから）プラマイをしっかり区別しないとならぬ
一度条件文の修正を忘れてこうなった
これはこれでいとをかし
脛節にかかる縦軸方向の力も正弦波みたいになった
（図にはLiftと書いてあるが、流体力以外も含むので厳密にはLiftではない）
修正後
いい感じ
またぶっ飛んで欲しくないので姿勢拘束と地面の抵抗は0にしてある
左上に経過時間を表示してみたが、値がキリ悪すぎる
1000倍してint型にしてミリ秒としようかな
さて、
coilSpringDamperCOMを読んでると、
// Accumulate the force for the restrained body
if( bodyIndex_ != 0 ) fx[bodyIndex_] += spatialVector( Rr & moment, Zero /*force*/ );
if( referenceID != 0 ) fx[referenceID] -= spatialVector( Rr & moment, Zero /*force*/ );
bodyにかかるモーメントに対してreactionには逆向きにモーメントがかかる
並進方向の力は必ずゼロ
下図のように腿節と脛節の角度差が初期角度差から変化すると、コイルばねは両側に戻そうとする力がかかるため、腿節と脛節には反対のモーメントがかかる
角度差なので、例えば腿節が回ると、ばねの自然長（というより自然角）は変化する
解析上ではもちろん角速度差に乗算されるダンパ係数もある
pythonコードは具体的にこんな感じ
### coilSpringDamper(tibia) ###
t_tibiaR = -ctrl.coilSpring * (ID_tibiaR.angPos-(ID_femurR.angPos-ID_femurR.initAng)) - ctrl.coilDamper * (ID_tibiaR.angVel-ID_femurR.angVel)
t_tibiaL = -ctrl.coilSpring * (ID_tibiaL.angPos-(ID_femurL.angPos-ID_femurL.initAng)) - ctrl.coilDamper * (ID_tibiaL.angVel-ID_femurL.angVel)
先も述べた通り、bodyとreactionで反対のモーメントがかかるので、
正負に気を付けて書き加えればよい
とりあえず、こんな値で動かしてみると
coilSpring = 5e-4
coilDamper = 1e-7
固すぎだがとりあえず動いた
で、バネダンパの調整がなんかうまくいかんなーと思ってたら、モーメントの計算がおかしいことに気づいた
# 回転
ID_body.angAcc = (-t_femurR_body - t_femurL_body) / ID_body.inertia
ID_femurR.angAcc = (t_femurR + t_femurR_body + t_femurR_tibiaR - t_tibiaR) / ID_femurR.inertia
ID_femurL.angAcc = (t_femurL + t_femurL_body + t_femurL_tibiaL - t_tibiaL) / ID_femurL.inertia
ID_tibiaR.angAcc = (t_tibiaR - t_femurR_tibiaR + t_tibiaR_ground) / ID_tibiaR.inertia
ID_tibiaL.angAcc = (t_tibiaL - t_femurL_tibiaL + t_tibiaL_ground) / ID_tibiaL.inertia
問題は並進ばねによるモーメント計算であり、
モーメントは剛体重心からの距離と、ある点にかかる力ベクトルとの外積で求まるため、
それぞれ長さの違う剛体同士で正負対称の値にはならない
言い訳としては、並進の計算をコピペして書いたためである（orz）
直した
# 回転
ID_body.angAcc = (t_body_femurR + t_body_femurL) / ID_body.inertia
ID_femurR.angAcc = (t_femurR + t_femurR_body + t_femurR_tibiaR - t_tibiaR) / ID_femurR.inertia
ID_femurL.angAcc = (t_femurL + t_femurL_body + t_femurL_tibiaL - t_tibiaL) / ID_femurL.inertia
ID_tibiaR.angAcc = (t_tibiaR_femurR + t_tibiaRsR_ground + t_tibiaRsL_ground + t_tibiaR) / ID_tibiaR.inertia
ID_tibiaL.angAcc = (t_tibiaL_femurL + t_tibiaLsR_ground + t_tibiaLsL_ground + t_tibiaL) / ID_tibiaL.inertia
それぞれの剛体が受けるモーメントで足りなかったものを新しく定義した
腿節脛節間の姿勢拘束をめっちゃ固くしたもの
ほぼ1関節モデル
姿勢拘束を調節したもの
それぞれの跳躍高度
離水時の速度が異なるため、位置の差は拡大していく
実際には、一度メニスカスを破壊すると肢は水没して上昇しようとする力に対して粘性がかかるのでさらに差は大きくなると考えられる
OpenFOAM（というかVOF法）では、低い計算コストで高精度な界面を表現するのがどうしても難しい。そのため、数値解析では実験や解析解よりも浅いところでメニスカス破壊が起きる。pythonで計算した通りに動いたとしても、OpenFOAMではメニスカス破壊が起きる可能性が高い