歯形係数の計算プログラムをExcelVBAで作る

経緯

「CADで危険断面位置と歯形係数を求める」では、実際の歯の形から、必要な寸法を計測して歯形係数を求めました。設計実務では、このような図形処理ではなく数値処理で求めるのですが、その内容を紹介します。

利用する歯形係数計算式

計算式の記載ある資料

金属歯車用の強度計算式のうち、歯形係数計算式が記載されているのは、国内に限ると以下の2つが挙げられます。

  • 日本歯車工業会JGMA401
  • 日本歯車工業会JGMA6101

プラスチック歯車では、次の資料に歯形係数計算式が載っています。

  • JIS B1759-2019 プラスチック円筒歯車の曲げ強さ評価方法

金属歯車用なので、ふさわしいのはJGMA401か、JGMA6101です。自分もこちらを使っているのですが、資料が市販品のため、プログラム部を公開するのは問題がありそうです。

いっぽう、プラスチック用JIS B 1759は、ネットで閲覧可能なので、こちらを利用することにします。
kikakurui.com

また、以下のサイトで詳細に解析されていますから、中身を知りたい方はお勧めです。
opeo.jp

JGMAとJISB1759の歯形係数の違い

両者は、「かみ合い率」の補正項の入れ方がことなります。これについて若干の説明をしましょう。
「かみ合い率」は平歯車では1以上2未満のような値を取ります。かみ合い長さと法線ピッチの比ですから、1歯だけを見ると、かみ合い率1.6を例にとると、1かみ合い長さが法線ピッチの1.6倍です。

図1.かみ合い率1.6の意味

次に先行する歯も同様に1かみ合い長さで1.6法線ピッチあり、1ピッチ先行するので0.6ピッチ重なります。

今度は後続する歯を見ると1ピッチ遅れるので、現行歯とは0.6ピッチ重なります。
つまりかみ合い率1.6の歯は、先行歯と0.6ピッチ、後続歯と0.6ピッチ重なるので、どことも重ならないのは1.6-0.6-0.6=0.4ピッチです。

言い換えると、かみ合い率1.6の歯車は、前後0.6ずつが2歯でかみ合い、残り0.4が1歯でかみ合います。

それでは、強度で問題になるのはどう言う状況かというと
歯先に荷重を受けると、モーメントが最大になるが
歯先は前後の2歯で分担するので、それぞれの歯の先端荷重は1/2になる
という理由で歯先荷重は最悪負荷条件ではないのです。

最悪荷重条件となるのは、1歯かみ合い中でもっともモーメントが大きくなる点、すなわち1歯かみ合いの最外側点です。

図2.歯先荷重点と最悪荷重点

この認識は一般的なのですが、計算への反映の仕方が、歯車規格毎に異なります。

計算式 歯形係数 補正項
JGMA401,6101 歯先荷重 あり
JIS B 1759 一歯最悪荷重 なし

つまり、JIS B 1957は歯形係数を最悪荷重点で計算し、JGMAは歯形係数は歯先で計算してから、「かみ合い率係数」という補正項を加算します。

どっちを使うか

このブログで何回か報告している歯形係数は「歯先荷重」のほうです。
そこで、計算式としてはJIS B 1759を使うが、作用点は歯先荷重として修正することにしました。

図.荷重点の違い

使用ツール

誰もがお手軽に計算できるように、Excel VBAを使いました。みなさんが自分でインターフェースを工夫されることを期待して、ここでは素のVBAだけを紹介します。

Excel「開発」「VisualBasic」で「モジュール」を挿入し、次のVBAを貼り付けて、singleCaseMain関数またはmultiCaseMain関数を実行させます
ツールバーに「開発」がない場合は、「オプション」「リボンのユーザー設定」で「開発」をリボンに表示させます。

JIS B 1759の計算式

以下にJIS B 1759に記載されている計算式をVBAに置きなおしたものです。
1行目のE式の右辺第3項が「-1(-sin(an))」というのはJISのミスプリントで、正しくは「-(1-sin(an))」です。

'JIS B 1759オリジナル

    E = pi / 4 - hfp * Tan(an) - 1(- Sin(an)) * ro / Cos(an) ---誤
    G = ro / mn - hfp / mn + x
    H = 2 / zn * (pi / 2 - E / mn) - pi / 3
    
 '次式をθについて解く
    th = 2 * G / zn * Tan(θ) - H
    
    
    '一歯かみ合いの外側径入力値
    den =()
    '一歯かみ合いの外側入力荷重角度
    afen =()
    '一歯かみ合いの外側入力値
    Ye =()
    '危険断面幅
    sfn = (zn * Sin(pi / 3 - th) + Sqr(3) * (G / Cos(th) - ro / mn)) * mn
    '作用点までの腕の長さ
     hfe = 1 / 2 * ((Cos(Ye) - Sin(Ye) * Tan(afen)) * den / mn 
                    - zn * Cos(pi / 3 - th) - (G / Cos(th) - ro / mn))
    
    '#歯形係数
    yf = (6 * hfe/mn * Cos(αFen)) / ((sfn/mn * sfn/mn) * Cos(wt))

作成したVBAプログラム

単発計算部 sub singleCaseMain()

1セットの諸元問い合わせに対して、その歯形係数を計算して表示します。
ユーザはエディタで次の諸元を変更して、実行します。

mn = 1 '歯直角モジュール
zn = 20 '歯数
an = 20 '圧力角
x = 0# '転位係数

hk = 1 '歯末のたけ(mnで規格化)
hfp = 1.25 '歯元のたけ(mnで規格化)
ro = 0.375 'ラック歯先の丸み(mnで規

「CADで危険断面位置と断面係数を求める」で、実測値と理論値の比較をしてきましが、ここでは理論値同士の比較をやってみましょう。

上記諸元結果
項目 今回理論値 前回理論値
危険断面厚さ 1.9440374 1.9438910
作用点までの長さ 1.9071318 1.9071201
歯形係数 2.804039 2.804444

差は0.014%なので、まったく問題ありません。

複数計算部 sub multiCaseMain()

よくある歯形係数グラフ作成用です。
an,hk,hfp,roを固定して、転位係数x16水準、歯数zn31水準をすべて計算して、結果を表として「sheet1」に出力します。
表の1行目は歯数ですが、2行目は歯数の逆数になっています。これは、歯形係数グラフの横軸を歯数の逆数にしている例が多いため、それに合わせたものです。

グラフは、転位係数をパラメータとした歯数・歯形係数のxy図で作成し、歯数の逆数を横軸に取って、縦軸の範囲を1.8~3.8に設定します。グラフの横軸目盛を反転させると曲線部分は完成ですが、横軸目盛ラベルは歯数の逆数になっています。ここのラベルとしては元の歯数を使うべきですが、Excelでは目盛の差し替えができませんでした。

それでもできるだけ合わせこんで、以前作成した歯形係数グラフと比較したのが下図です。グラフを重ねるとぴったり一致するので、全領域計算は合っているようです。

図.歯形係数グラフの比較

歯形係数計算部 private function toothFormFactor()

オリジナルに対して次を変更しました。

  • ro,hfpをmnで規格化
  • 一歯かみ合いの外側荷重入力に対して、歯先荷重入力に変更
  • θの収束計算にニュートンラフソン法の導入

求めたいのは、次式を満足するθなのですが、右辺にtanθが含まれるため、解析的には求めらえません。

θ=2 * G / zn * tanθ - H

こういう場合に使うのは二分法とかニュートンラフソン法です。なかでもニュートンラフソン法は収束が早く、今回の例でも数回の繰り返しで収束します。

  • 応力集中係数は割愛

「CADで危険断面位置と断面係数を求める」で求めたレベルの歯形係数に留めました。

ソースコード全文

Option Explicit
'deg-rad変換
Private Function radians(v)
    radians = v * WorksheetFunction.pi() / 180
End Function
'rad-deg変換
Private Function degrees(v)
    degrees = v / WorksheetFunction.pi() * 180
End Function
'inv関数
Private Function inv(v)
    inv = Tan(v) - v
End Function
'単発計算用
Sub singleCaseMain()
    Dim mn#, zn#, an#, ro#
    Dim hfp#, x#, hk#, yf#
    
    '諸元入力
    mn = 1              '歯直角モジュール
    zn = 20             '歯数
    an = 20             '圧力角
    x = 0#              '転位係数
    
    hk = 1              '歯末のたけ(mnで規格化)
    hfp = 1.25          '歯元のたけ(mnで規格化)
    ro = 0.375          'ラック歯先の丸み(mnで規
    
    '歯形係数計算
    yf = toothFormFactor(mn, zn, an, x, hk, hfp, ro)
    
    Debug.Print yf
End Sub

'歯数-転位係数グラフのための作表

Sub mutiCaseMain()
    Dim z_array
    Dim x_array
    Dim sht
    
    '"sheet1"シートがないとエラー
    Set sht = Sheets("sheet1")
    sht.Range("A1:AZ100").ClearContents
    
    '転位係数配列
    x_array = Array(-0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, _
        0.5, 0.6, 0.7, 0.8, 0.9, 1#)
    '歯数配列
    z_array = Array(6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, _
        30, 40, 50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 160, 170, _
        180, 190, 200, 400, 1000)
    
    '歯数と転位係数の組み合わせ計算
    Dim x_idx%, z_idx%, x#, z#, yf#
    For x_idx = LBound(x_array) To UBound(x_array)
        '転位係数見出し
        sht.Cells(x_idx + 3, 1) = x_array(x_idx)
        
        For z_idx = LBound(z_array) To UBound(z_array)
            If x_idx = 1 Then
                ' 歯数見出し
                sht.Cells(1, z_idx + 2) = z_array(z_idx)
                sht.Cells(2, z_idx + 2) = 1 / z_array(z_idx)
            End If
            
            x = x_array(x_idx)
            z = z_array(z_idx)
            
            ' 歯形係数計算
            yf = toothFormFactor(1, z, 20, x, 1#, 1.25, 0.375)
            
            'シート記載
            sht.Cells(x_idx + 3, z_idx + 2) = yf
         
        Next
     Next
End Sub

'歯形係数計算
Private Function toothFormFactor(mn#, zn#, an#, x#, hk#, hfp#, ro#)
    Dim E#, H#, G#
    Dim th#, thNew#
    Dim F#, df#
    Dim pi#
    
    pi = WorksheetFunction.pi()
    an = radians(an)
    
    ' ラック歯先の丸み中心x位置
    E = (pi / 4 - hfp * Tan(an) - (1 - Sin(an)) * ro / Cos(an)) * mn
    
    ' ラック歯先の丸み中心y位置
    G = (ro - hfp + x)
    
    H = (2 / zn * (pi / 2 - E / mn) - pi / 3)
 
    'NewtonRaphson法で30°接線を求める
    Dim N!
    
    th = pi / 6     '角度初期値
    thNew = 0       '角度更新値
    N = 0
    Do While (N < 10)
        '目的関数
        F = 2 * G / zn * Tan(th) - H - th
        
        '1次導関数
        df = 2 * G / zn * (Tan(th) ^ 2 + 1) - 1
        
        'NewtonRaphson法
        thNew = th - F / df
        
        '収束したらループ抜ける
        If Abs(thNew - th) < 0.000001 Then
            Exit Do
            
        '収束しない場合は値を更新してループ
        Else
            th = thNew
        
        End If
        N = N + 1

    Loop
    '30°接線となる角度th
    th = thNew
    
    '以降は
    'JIS B 1759の歯形係数は、一歯かみ合い領域の最外側で荷重負荷を考える
    'JGMAの歯形係数は、二歯かみ合いの歯先荷重負荷で、かみ合い率は補正項で考慮する
    
    'このプログラムは、30°接線位置の算出にJIS B 1759記載の計算方法を利用したが、
    '荷重負荷は、歯先方式とする。

    Dim AK#, Delta#, anf#, da#
    
    ' 歯先圧力角
    AK = WorksheetFunction.Acos(zn * Cos(an) / (zn + 2 * (hk + x)))
    ' 歯先歯厚角の1/2
    Delta = 1 / zn * (pi / 2 + 2 * x * Tan(an)) + inv(an) - inv(AK)
    ' 歯直角歯先荷重角
    anf = AK - Delta
    ' 歯先径
    da = (zn + 2 * (hk + x)) * mn
    
    Dim sfn#, hfe#, Ye#, den#, afen#
    ' den 相当平歯車の一歯かみ合い領域の歯先側の点を通る円の直径(mm)
    'γe :相当平歯車の一歯かみ合い領域の歯先側の点の角度(値表示は°,計算はrad)

    '一歯かみ合いの外側径入力値に、歯先径を代入
    den = da
    
    '一歯かみ合いの外側入力荷重角度に、歯先への荷重角度代入
    afen = anf
    
    '一歯かみ合いの外側入力値に、歯先歯厚の半角を代入
    Ye = Delta
    
    '危険断面幅:規格化したro使用
    sfn = (zn * Sin(pi / 3 - th) + Sqr(3) * (G / Cos(th) - ro))
    
    '作用点までの腕の長さ:規格化したro使用
    hfe = 1 / 2 * ((Cos(Ye) - Sin(Ye) * Tan(afen)) * den / mn - _
            zn * Cos(pi / 3 - th) - (G / Cos(th) - ro))
    
    'Debug.Print "危険断面幅"; , sfn * mn, "作用点までの腕の長さ:", hfe * mn; ""
    
    '#歯形係数
    Dim yf#
    yf = (6 * hfe * Cos(anf)) / ((sfn * sfn) * Cos(an))
    
    'Debug.Print "歯形係数", yf
    
    toothFormFactor = yf
    
End Function