MV-9001：シンプルな振子のチュートリアル

1. このコマンドを実行します：

   In [3]: units= Units () 

2. mmksでモデル化したいので、単位のlengthプロパティを変更し、"MILLIMETERS"の値を割り当てます。簡単に直接割り当てることができます。

   In [4]: units.length="MILLIMETER"

3. 同様に、Accgravクラスを使って重力加速度場を作ることができます。

   In [5]: gravity=Accgrav(kgrav=-9810)

   Accgravの作成結果は、gravityという変数に格納されます。それが等号の左に出てくる名称です。これは現在のネームスペース内のPython変数であり、スコープ内のどこでも使用できます。

msolveのクラス名ではCapWordsの規約を使用し、Pythonのヘルプシステムを使用して各msolveクラス、その使用方法、プロパティについて学ぶことができます。例として以下の手順を実行してください。

4. Accgravクラスと対話するには、クラス名を入力引数として、help ( )関数を呼び出します。

   In [6]: help(Accgrav)
   Help on class Accgrav in module msolve.Model:
   class Accgrav(Force)
    |  Accgrav defines the acceleration due to gravity along the global X, Y, and Z directions.
    |
    |  Method resolution order:
    |      Accgrav
    |      Force
    |      msolve.Base.SolverObject
    |      __builtin__.object
    |
    |  Methods defined here:
    |
    |  sendToSolver(self)
    |      Send the object to the solver
    |
    |  setSolverValue(self, name, value)
    |      Set the property 'name' to 'value' between simulations
    |
    |  ----------------------------------------------------------------------
    |  Data descriptors defined here:
    |
    |  igrav
    |      Gravity in the global X direction
    |      Type=Double, Optional, Default=0
    |      Modifiable during simulation
    |
    |  jgrav
    |      Gravity in the global Y direction
    |      Type=Double, Optional, Default=0
    |      Modifiable during simulation
    |
    |  kgrav
    |      Gravity in the global Z direction
    |      Type=Double, Optional, Default=0
    |      Modifiable during simulation
   

   ヘルプ関数が呼び出されると、多くの情報が表示されます。情報の中で最も重要なのは、クラスとそのプロパティの記述です。例えば、データ記述子は以下のコードで定義されています：

   igrav
       Gravity in the global X direction
       Type=Double, Optional, Default=0
       Modifiable during simulation

   msolveでは、クラスのインスタンスを作成すると必ずオブジェクトが返されます。これにより、オブジェクトの名称、この場合は重力を使用して参照する機能が得られます。 重力とその特性はモデル内のどこでも使用できます。例えば、モデル構築の段階、もしくはMotionSolveシミュレーションの間でプロパティを変更することができます。このアプローチの利点については、本チュートリアル内で後程説明します。

   一般的なルールとして、作成されたオブジェクトをローカル変数を使って保存するのは有益です。この時点で、モデルの作成に進むことができます。

5. 地面のパートを作成し、Partクラスを使用して、地面のBoolean属性を'True’に指定します。

   In [7]: ground=Part(ground=True) 

1. 以下のコマンドで、 2つの軸のそれぞれの方向付けに用いられる空間内のポイントを与えると、3番目は直交右手座標系として自動的に計算されます。

   In [9]: ground_mar=Marker (body=ground, qp=p0, zp=[0,100,0], xp=[100,0,0]) 

2. これで、振子ボディを作成することができます。質量が2kgで、対角項のみに慣性特性のある剛体パートです。

   In [10]: part=Part(mass=2.0, ip=[1e3,1e3,1e3]) 

3. オプション: コードセルにpartとTABキーを入力すると、partに適用されるすべての機能とプロパティがドロップダウンリストとして表示されます。

1. validate関数を呼び出してパートを検証します：

   In [11]: part.validate() 
   Part/2
   ERROR:: Mass is specified but cm is not specified.
   ERROR:: There are no markers on this part.
   Out[11]:
   False

   エラーメッセージは、パートの重心（cm）が欠落していることを告げています。
2. この問題を修正するには、(100, 0, 0)の位置にパートの重心マーカーを作成します。このポイントは、振子パートの位置を表しています。

   In [12]: p1 = Point(100,0,0)
   # create the marker as the part.cm 
   # note that the '.' operator allows us to define a the 'cm' property of the part object. 
   # the marker will be positioned at (100,0,0), the X axis is pointing towards the point (200, 0,0) 
   # and the z axis points to (200, 100, 0)
    
   part.cm = Marker(body=part, qp=p1, zp=[200,100,0], xp=[200,0,0]) 
   

1. 以下のコマンドを実行します：

   In [14]: joint = Joint(type="REVOLUTE", i=ground_mar, j=Marker(body=part, qp=p0, xp=[100,0,0], zp=[0,100,0]))
   force_req = Request(type="FORCE", i=joint.i, j=joint.j, rm=joint.i)
   angle_req = Request(f1="RTOD*AX({I}, {J})".format(I=joint.i.id, J=joint.j.id),                    
                       f2="RTOD*AY({I}, {J})".format(I=joint.i.id, J=joint.j.id),
                       f3="RTOD*AZ({I}, {J})".format(I=joint.i.id, J=joint.j.id)) 

2. Sforceを呼び出して、ボールが動くときの接触力を計算します。入力関数はMotionSolve式である必要があり、マーカーjのz軸は力がかかる方向を定義します。

   In [15]: input_function =
   "IMPACT(DZ({I}, {J}, {J}),VZ({I}, {J}, {J}),{G},{K},{E},{C},{D})"
   .
   format( I=marker_i.id,J=marker_j.id, K=5e9, G=0.00, E=1.0, C=0.8e5, D=0.0 )
   contact_force = Sforce(i = marker_i.id, j = marker_j.id,
   type=
   "TRANSLATION"
   , function = input_function) 
   

1. 以下のコマンドを実行します：

   In [16]: run=model.simulate (type="TRANSIENT", returnResults=True, end=2, dtout=.01)

2. 直接代入で値を変更することで、パートの質量プロパティを簡単に変更することができます。

   In [17]: part.mass= 12

3. この時点で、前回の解析の最後から別のシミュレーションを実行することができます。これは単に、新しいシミュレートコマンドを発行し、終了時間を4秒に指定するだけで行えます。なお、このシミュレーションは新しいパートの質量で行われているため、その影響をジョイントの反力で見ることができるはずです。

   In [18]: run=model.simulate (type="TRANSIENT", returnResults=True, end=4, dtout=.01)

4. 次のコマンドで、実行コンテナからforce_reqを抽出しています。時間変化する結果を抽出するのに便利なデータ構造を持つことが目的です。

   In [19]: force=run.getObject (force_req)

5. 力のオブジェクトを使用して、縦方向の反力Fxをプロットすることができます。そのために、IPythonノートブックで利用可能なmatplotlibモジュールを使用しています。プロットはインラインで表示されます。時間＝2秒での質量変化が反力の大きさに影響を与え、350Nのピークを発生させていることがわかります。

   In [20]: %matplotlib inline
   from matplotlibimport
   pyplot pyplot.plot (force.times, force.getComponent (2), label=force.labels[2])

   
   
   図 2. Out[20]: [<matplotlib.lines.Line2D at 0x7142860>]

6. 同様に、回転角をプロットすることができます。角変位リクエストでは一般的なランタイム式が使われていたため、AY（角リクエストのマーカjのY軸周りの角度）に対応する正しい成分を取得する必要があります。

   In [21]: angle = run.getObject(angle_req)
   pyplot.plot (angle.times, angle.getComponent(2), label="angle [deg]")

   
   
   図 3. Out[21]: [<matplotlib.lines.Line2D at 0x72c8fd0>]

   角変位は単純な調和関数であり、接合部に粘性摩擦などの散逸力がないため、振子は180度に達します。振動の周波数はこれらの式に支配されていることにご留意ください。(1) $$T=2\pi \sqrt{\frac{L}{g}}$$ (2) $$f=\frac{1}{T}$$

   上のプロットから分かるように、質量は約1.5Hzの周波数で揺れています。以下のgenerateOutput()関数は、HyperViewとHyperView Playerで可視化するためのMotionSolve *.pltファイルとH3Dアニメーションファイルを生成します。

   In [22]: model.generateOutput()