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

2020/03/10 categories:Python| tags:Python|

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

NumPyのmeshgridで総当たり計算

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

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

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

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

実行結果

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

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

print( c.flatten() )
print( d.flatten() )

実行結果

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

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

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

実行結果

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

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

クラスから関数へ

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

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

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

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

実行結果

[11629.62962963  9303.7037037  11629.62962963  9303.7037037
 12717.         10173.6        12717.         10173.6       ]

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

NumPyでばねの計算

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

d = np.arange(0.2, 0.5, 0.1) # 線径
G = np.array([78500.0]) # 横弾性係数
Di = np.arange(5.0, 8.0, 1.0) # 内径
Na = np.arange(5.0, 6.0, 1.0) # 有効巻数
Ns = np.arange(1.0, 2.0, 1.0) # 座の巻き数
Ho = np.arange(10.0, 12.0, 1.0) # 自由長さ
H1 = np.arange(8.0, 9.0, 1.0) # 最小荷重時長さ
H2 = np.arange(4.0, 5.0, 1.0) # 使用時長さ
no_g = np.array([1.0]) # 研削
(d,G,Di,Na,Ns,Ho,H1,H2,no_g) = np.meshgrid(d,G,Di,Na,Ns,Ho,H1,H2,no_g)
(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())
_tuple = coil_spring_numpy(d,G,Di,Na,Ns,Ho,H1,H2,no_g)
datas = np.array(_tuple).T

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

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

ソースコード

## -*- coding: utf-8 -*-
import math
import numpy as np
def coil_spring_numpy(d,G,Di,Na,Ns,Ho,H1,H2,no_grinding):
    Do = Di + d*2.0 # 外形
    Nt = Na + Ns*2.0 # 総巻数
    D = Di + d # コイル平均径
    ar = Ho / D # 縦横比
    c = D / d # ばね指数
    k = G * (d**4.0) / (8.0 * Na * (D**3.0)) # ばね定数
    p = (Ho - (Nt - Na + no_grinding) * d) / Na # ピッチ
    K = (4.0 * c - 1.0) / (4.0 * c - 4.0) + 0.615 / c # 応力修正係数
    Hs = (Nt + 1.0) * d # 密着長さ
    Fs = k * (Ho - Hs) # 密着荷重
    Rs = (Ho - Hs) / (Ho - Hs) * 100.0 # 動作範囲
    TSs = 8.0 * D * Fs / (math.pi * (d**3.0)) # ねじり応力
    TCSs = K * TSs # 密着ねじり修正応力
    F1 = k * (Ho - H1) # 最小荷重
    R1 = (Ho - H1) / (Ho - Hs) * 100.0 # 動作範囲
    TS1 = 8.0 * D * F1 / (math.pi * (d**3.0)) # ねじり応力
    TCS1 = K * TS1 # 最小ねじり修正応力
    F2 = k * (Ho - H2) # 使用時荷重
    R2 = (Ho - H2) / (Ho - Hs) * 100.0 # 動作範囲
    TS2 = 8.0 * D * F2 / (math.pi * (d**3.0)) # ねじり応力
    TCS2 = K * TS2 # 使用時ねじり修正応力
    H3 = 0.2*Ho + 0.8*Hs # 最大荷重時長さ
    F3 = k * (Ho - H3) # 最大荷重
    R3 = (Ho - H3) / (Ho - Hs) * 100.0 # 動作範囲
    TS3 = 8.0 * D * F3 / (math.pi * (d**3.0)) # ねじり応力
    TCS3 = K * TS3 # 最大ねじり修正応力
    return (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)
if __name__ == '__main__':
    d = np.arange(0.2, 0.5, 0.1) # 線径
    G = np.array([78500.0]) # 横弾性係数
    Di = np.arange(5.0, 8.0, 1.0) # 内径
    Na = np.arange(5.0, 6.0, 1.0) # 有効巻数
    Ns = np.arange(1.0, 2.0, 1.0) # 座の巻き数
    Ho = np.arange(10.0, 12.0, 1.0) # 自由長さ
    H1 = np.arange(8.0, 9.0, 1.0) # 最小荷重時長さ
    H2 = np.arange(4.0, 5.0, 1.0) # 使用時長さ
    no_g = np.array([1.0]) # 研削
    
    (d,G,Di,Na,Ns,Ho,H1,H2,no_g) = np.meshgrid(d,G,Di,Na,Ns,Ho,H1,H2,no_g)
    
    (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())
    _tuple = coil_spring_numpy(d,G,Di,Na,Ns,Ho,H1,H2,no_g)
    datas = np.array(_tuple).T
    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"]
    dicts = [ { key:round(d,6) for key, d in zip(keys, data) } for data in datas ]
    print(dicts)