Pythonで圧縮コイルバネの計算をしてみる NumPy版

Share on:

以前ばねの設計用ツールを作るためにSpringというクラスを作成して、ばね関連の計算を行うプログラムを作成しました。しかし、そのクラスを使ってツールを作ると実行速度がかなり遅くなってしまったため、NumPyを用いて計算の高速化をすることにしました。

NumPyのmeshgridで総当たり計算

圧縮コイルバネの計算をするにあたって下記のパラメータを幅を持たせて計算しようと思います。

  • 線径
  • 横弾性係数
  • 内径
  • 有効巻き数
  • 座の巻き数
  • 自由長さ
  • 最小荷重時長さ
  • 使用時長さ
  • 切削の有無

これらのパラメータを全通り計算しようとしたときにNumPyのmeshgridを使うといいみたいです。(全通り計算する必要があるかはとりあえず置いておきます。)例えば、下記のようにnp.arrayのa、bのmeshgridを求めてみます。

1import numpy as np
2a = np.array([1, 5, 10])
3b = np.array([2, 6,  3])
4c, d = np.meshgrid(a, b)
5print(a)
6print(b)
7print(c)
8print(d)

実行結果

1[ 1  5 10]
2[2 6 3]
3[[ 1  5 10]
4 [ 1  5 10]
5 [ 1  5 10]]
6[[2 2 2]
7 [6 6 6]
8 [3 3 3]]

これらc、dに対してNumPyのflattenを実行してあげると、

1print( c.flatten() )
2print( d.flatten() )

実行結果

1[ 1  5 10  1  5 10  1  5 10]
2[2 2 2 6 6 6 3 3 3]

このc.flatten、d.flattenをarrayとして行列の転置をすると

1print( np.array( (c.flatten(), d.flatten())).T )

実行結果

1[[ 1  2]
2 [ 5  2]
3 [10  2]
4 [ 1  6]
5 [ 5  6]
6 [10  6]
7 [ 1  3]
8 [ 5  3]
9 [10  3]]

というように、組み合わせの配列が得られるので、これを前述のパラメータに適用して全通り計算をしてみようと思います。

クラスから関数へ

NumPyで計算するために、ばねクラスを関数に書き換えました。例えば、ばね定数を求める計算は、

1k = G * (d**4.0) / (8.0 * Na * (D**3.0)) # ばね定数

というように単純な計算の式ですが、この変数をそのままNumPyのarrayに置き換えてしまえばそのまま計算できるので非常に楽です。例えば、Gを1種類、dを2種類、Naを2種類、Dを2種類で計算するとした場合、

1G  = np.array([78500.0, 78500.0, 78500.0, 78500.0, 78500.0, 78500.0, 78500.0, 78500.0])
2d  = np.array([2.0, 2.0, 2.0, 2.0, 3.0, 3.0, 3.0, 3.0])
3Na = np.array([4.0, 5.0, 4.0, 5.0, 4.0, 5.0, 4.0, 5.0])
4D  = np.array([1.5, 1.5, 1.5, 1.5, 2.5, 2.5, 2.5, 2.5])
5k  = G * (d**4.0) / (8.0 * Na * (D**3.0))
6print(k)

実行結果

1[11629.62962963  9303.7037037  11629.62962963  9303.7037037
2 12717.         10173.6        12717.         10173.6       ]

こんな感じで簡単にNumPyのarrayだけで、Pythonのforを使わずに計算できてしまい、計算速度がかなり速くなります。

NumPyでばねの計算

ばねの計算をcoil_spring_numpyというメソッドで計算するとすると

 1d = np.arange(0.2, 0.5, 0.1) # 線径
 2G = np.array([78500.0]) # 横弾性係数
 3Di = np.arange(5.0, 8.0, 1.0) # 内径
 4Na = np.arange(5.0, 6.0, 1.0) # 有効巻数
 5Ns = np.arange(1.0, 2.0, 1.0) # 座の巻き数
 6Ho = np.arange(10.0, 12.0, 1.0) # 自由長さ
 7H1 = np.arange(8.0, 9.0, 1.0) # 最小荷重時長さ
 8H2 = np.arange(4.0, 5.0, 1.0) # 使用時長さ
 9no_g = np.array([1.0]) # 研削
10(d,G,Di,Na,Ns,Ho,H1,H2,no_g) = np.meshgrid(d,G,Di,Na,Ns,Ho,H1,H2,no_g)
11(d,G,Di,Na,Ns, Ho,H1,H2,no_g) = (d.flatten(), G.flatten(), Di.flatten(), Na.flatten(), 
12    Ns.flatten(), Ho.flatten(), H1.flatten(), H2.flatten(), no_g.flatten())
13_tuple = coil_spring_numpy(d,G,Di,Na,Ns,Ho,H1,H2,no_g)
14datas = np.array(_tuple).T

というように、meshgridとflattenでばねのパラメータの組み合わせを作って、それをcoil_springに渡して、最後Tで転置したら結果が得られます。

ただし、ばねのパラメータの数を増やしまくったら実行時にPCがフリーズしたので配列を分割して実行したほうがいいみたいです。

ソースコード

 1# -*- coding: utf-8 -*-
 2import math
 3import numpy as np
 4def coil_spring_numpy(d,G,Di,Na,Ns,Ho,H1,H2,no_grinding):
 5    Do = Di + d*2.0 # 外形
 6    Nt = Na + Ns*2.0 # 総巻数
 7    D = Di + d # コイル平均径
 8    ar = Ho / D # 縦横比
 9    c = D / d # ばね指数
10    k = G * (d**4.0) / (8.0 * Na * (D**3.0)) # ばね定数
11    p = (Ho - (Nt - Na + no_grinding) * d) / Na # ピッチ
12    K = (4.0 * c - 1.0) / (4.0 * c - 4.0) + 0.615 / c # 応力修正係数
13    Hs = (Nt + 1.0) * d # 密着長さ
14    Fs = k * (Ho - Hs) # 密着荷重
15    Rs = (Ho - Hs) / (Ho - Hs) * 100.0 # 動作範囲
16    TSs = 8.0 * D * Fs / (math.pi * (d**3.0)) # ねじり応力
17    TCSs = K * TSs # 密着ねじり修正応力
18    F1 = k * (Ho - H1) # 最小荷重
19    R1 = (Ho - H1) / (Ho - Hs) * 100.0 # 動作範囲
20    TS1 = 8.0 * D * F1 / (math.pi * (d**3.0)) # ねじり応力
21    TCS1 = K * TS1 # 最小ねじり修正応力
22    F2 = k * (Ho - H2) # 使用時荷重
23    R2 = (Ho - H2) / (Ho - Hs) * 100.0 # 動作範囲
24    TS2 = 8.0 * D * F2 / (math.pi * (d**3.0)) # ねじり応力
25    TCS2 = K * TS2 # 使用時ねじり修正応力
26    H3 = 0.2*Ho + 0.8*Hs # 最大荷重時長さ
27    F3 = k * (Ho - H3) # 最大荷重
28    R3 = (Ho - H3) / (Ho - Hs) * 100.0 # 動作範囲
29    TS3 = 8.0 * D * F3 / (math.pi * (d**3.0)) # ねじり応力
30    TCS3 = K * TS3 # 最大ねじり修正応力
31    return (d,G,Di,Na,Ns,Ho,H1,H2,no_grinding,Do,Nt,D,ar,c,k,p,K,Hs,Fs,Rs,TSs,TCSs,
32        F1,R1,TS1,TCS1,F2,R2,TS2,TCS2,H3,F3,R3,TS3,TCS3)
33if __name__ == '__main__':
34    d = np.arange(0.2, 0.5, 0.1) # 線径
35    G = np.array([78500.0]) # 横弾性係数
36    Di = np.arange(5.0, 8.0, 1.0) # 内径
37    Na = np.arange(5.0, 6.0, 1.0) # 有効巻数
38    Ns = np.arange(1.0, 2.0, 1.0) # 座の巻き数
39    Ho = np.arange(10.0, 12.0, 1.0) # 自由長さ
40    H1 = np.arange(8.0, 9.0, 1.0) # 最小荷重時長さ
41    H2 = np.arange(4.0, 5.0, 1.0) # 使用時長さ
42    no_g = np.array([1.0]) # 研削
43    
44    (d,G,Di,Na,Ns,Ho,H1,H2,no_g) = np.meshgrid(d,G,Di,Na,Ns,Ho,H1,H2,no_g)
45    
46    (d,G,Di,Na,Ns, Ho,H1,H2,no_g) = (d.flatten(), G.flatten(), Di.flatten(), Na.flatten(), Ns.flatten(), Ho.flatten(), H1.flatten(), H2.flatten(), no_g.flatten())
47    _tuple = coil_spring_numpy(d,G,Di,Na,Ns,Ho,H1,H2,no_g)
48    datas = np.array(_tuple).T
49    keys = ["d","G","Di","Na","Ns","Ho","H1","H2","no_grinding","Do","Nt","D","ar","c","k","p","K","Hs","Fs","Rs","TSs","TCSs","F1","R1","TS1","TCS1","F2","R2","TS2","TCS2","H3","F3","R3","TS3","TCS3"]
50    dicts = [ { key:round(d,6) for key, d in zip(keys, data) } for data in datas ]
51    print(dicts)

関連記事