JP3910825B2 - Simulation method for performance prediction of products made of viscoelastic materials - Google Patents

Simulation method for performance prediction of products made of viscoelastic materials Download PDF

Info

Publication number
JP3910825B2
JP3910825B2 JP2001341117A JP2001341117A JP3910825B2 JP 3910825 B2 JP3910825 B2 JP 3910825B2 JP 2001341117 A JP2001341117 A JP 2001341117A JP 2001341117 A JP2001341117 A JP 2001341117A JP 3910825 B2 JP3910825 B2 JP 3910825B2
Authority
JP
Japan
Prior art keywords
strain
viscoelastic material
strain rate
viscoelastic
coordinate system
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
JP2001341117A
Other languages
Japanese (ja)
Other versions
JP2003139668A (en
Inventor
勝彦 植田
和佳 宮本
正貴 白石
Original Assignee
Sriスポーツ株式会社
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Sriスポーツ株式会社 filed Critical Sriスポーツ株式会社
Priority to JP2001341117A priority Critical patent/JP3910825B2/en
Priority to US10/481,754 priority patent/US7130748B2/en
Publication of JP2003139668A publication Critical patent/JP2003139668A/en
Application granted granted Critical
Publication of JP3910825B2 publication Critical patent/JP3910825B2/en
Anticipated expiration legal-status Critical
Expired - Fee Related legal-status Critical Current

Links

Images

Landscapes

  • Investigating Strength Of Materials By Application Of Mechanical Stress (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Description

【0001】
【発明の属する技術分野】
本発明は、粘弾性材料からなる製品の性能予測のためのシミュレーション方法に関し、詳しくは、粘弾性材料からなる製品の性能をシミュレーションにより精度良く予測するシミュレーション方法に関するものである。
【0002】
【従来の技術】
ゴムやエラストマーなどの高分子材料に代表される粘弾性材料は、タイヤ、スポーツ競技において使用される各種ボール、印刷機に用いられるロールなど各種製品に広く適用されている。
【0003】
上記のような粘弾性材料や金属材料等の各種工業製品について、試作の費用と時間の節約等のためにシミュレーションを用いた製品開発が多方面で行われている。例えば、ゴルフボールの反発性能を予測するために、実際の打撃試験のシミュレーション方法が提案されている。
【0004】
従来、上記のようなシミュレーションを行う場合、一般に材料の剛性、粘性などを測定する粘弾性スペクトルメーター、縦弾性係数(ヤング率)を測定する引張試験機などにより測定したボールの構成材料の物性値をシミュレーションへの入力値として用いている。特に、粘弾性スペクトルメーターでは、動的ひずみが与えられた試験片の物性値が測定されるため、粘弾性材料からなる製品のシミュレーションに有用である。
【0005】
【発明が解決しようとする課題】
しかしながら、上記した粘弾性スペクトルメーター、縦弾性係数を測定する引張試験機などによる測定では試験片に大きな変形量を与えることができず、測定時に、粘弾性材料からなる試験片にかかる最大ひずみ速度は、0.001/sから1.0/s程度と小さく、また、最大圧縮ひずみも0.0001から0.02程度と小さな値となる。
【0006】
一方、粘弾性材料からなる製品は、実使用時に外力の影響等で高速かつ大きな変形を起こす場合がある。例えば、ゴルフボールを構成する材料では、実際の打撃時に、最大ひずみ速度は500/s〜5000/s程度であり、最大圧縮ひずみは主に、0.05〜0.50程度であり、非常に大きな値となる。
【0007】
このように、粘弾性スペクトルメーター、縦弾性係数を測定する引張試験機等による測定では、上記のような実使用条件と同様の高速かつ大きな変形条件下での粘弾性材料の物性値を測定することはできず、実使用時と、最大ひずみ速度及び最大圧縮ひずみが大きく異なる。このため、粘弾性スペクトルメーター、引張試験機等により得られた材料物性値を入力する従来のシミュレーションでは、粘弾性材料物性を考慮した正確なシミュレーションを行うことができないという問題がある。
【0008】
即ち、衝撃荷重を受ける粘弾性材料の変形挙動は、静荷重を受ける場合とは異なり、変形量、あるいは変形速度によって著しく大きな影響を受けることが知られている。特にゴムやエラストマーなどに代表される高分子材料では、衝撃荷重により1万分あるいは1千分の数秒という非常に高速の変形を引き起こし、且つ変形量も数十%と非常に大きくなる。このような高速大変形の挙動を伴う粘弾性材料は多く、製品開発を効率的に行う上では精度の高いシミュレーションが要求されている。具体的には、ゴルフボール等のように使用時に衝撃を伴う製品では高速且つ大変形条件下での動的挙動および材料の特性が製品の性能を左右するため、実使用条件下での精度良いシミュレーションが製品開発に必要不可欠である。
【0009】
また、粘弾性材料には、衝撃荷重等の外力が加えられた場合、損失係数等の物性値がひずみ、ひずみ速度の大きさにより変化する材料がある。即ち、粘弾性材料は、変形の速度、大きさが広範囲であるため、変形速度や大きさによって、物性値が一次直線的でなく、強い非線形的な変化を伴う性質がある。具体的には、外力により材料が変形し、ひずみが大きくなるに伴い、S−S(ひずみ−応力)曲線のループ面積が大きくなり、損失係数等の物性値が、材料の変形状態(変形の速度、大きさ)等により変化し非線形性を示す。粘弾性材料によっては、強い非線形性を持つ材料も多く、このような粘弾性材料を用いた製品についても精度の高いシミュレーションが要求されている。
【0010】
しかしながら、粘弾性材料の非線形性を示す物性、例えば、損失係数が材料の変形の速度、大きさによって強く非線形的に変化する現象を精度良く表現できる手法がなかった。従来、粘弾性材料で構成されるゴルフボール等の製品については、粘弾性物性値がほとんど変化しないものとしてシミュレーションを行っていた。その結果、粘弾性材料からなる製品の実使用時の性能をシミュレーションにより正確に予測することができないという問題があり、実際には、試作品を製作することにより製品の性能を評価せざるを得ないという問題がある。
【0011】
本発明は上記した問題に鑑みてなされたものであり、材料物性が非線形性を示す粘弾性材料からなる製品の実使用条件下での性能を、精度良くシミュレーションにより予測することを課題としている。
【0012】
【課題を解決するための手段】
上記課題を解決するため、本発明は、粘弾性材料からなる製品の実使用状態を想定した測定条件下で、上記粘弾性材料に生じるひずみ、ひずみ速度、応力の時々刻々の値を測定し、
上記ひずみ、ひずみ速度、応力の時刻歴データと、上記粘弾性材料の粘性を考慮した粘弾性モデルとから、ひずみの負荷状態時とひずみの除荷あるいは復元状態時とで場合分けして上記粘弾性材料の粘性抵抗の時刻歴データを導出し、
上記粘弾性材料からなる製品を、解析対象の製品モデルとして設定し、該製品モデルを多数の要素に分割し、上記ひずみ、ひずみ速度、粘性抵抗の関係を入力し、上記ひずみの負荷状態時とひずみの除荷あるいは復元状態時とでの粘性抵抗の違いを考慮して有限要素法により解析を行い、上記粘弾性材料からなる製品モデルの性能を予測する際に、
上記要素毎に生じる全体座標系のひずみ及びひずみ速度を偏差成分と体積成分に分解し、上記偏差成分のひずみ及びひずみ速度を、全体座標系から主ひずみ座標系及び主ひずみ速度座標系に変換し、
上記変換した偏差主ひずみ及び偏差主ひずみ速度を用い、主ひずみ座標系及び主ひずみ速度座標系の各座標軸毎に、ひずみが同じ値でもひずみの負荷状態時とひずみの除荷あるいは復元状態時とで異なった粘性抵抗を決定していることを特徴とする粘弾性材料からなる製品の性能予測のためのシミュレーション方法を提供している。
【0013】
このように、実使用条件下で測定されたひずみ、ひずみ速度を体積成分と偏差成分に分解し、主ひずみ座標系及び主ひずみ速度座標系に変換した偏差成分の主ひずみ及び主ひずみ速度を用いて主ひずみ座標系及び主ひずみ速度座標系の各軸毎に、ひずみが同じ値の時でも、ひずみの負荷状態時とひずみの除荷あるいは復元状態時とで異なった粘性抵抗を決定し、粘性抵抗を要素毎に各軸方向3成分について可変にしている。このため、材料の変形状態が負荷状態であるのか、除荷状態であるのか、実際の状況に対応した精度良い解析を行うことができると共に、3次元方向について正確に解析を行うことができ、シミュレーションの性能予測精度を向上することができる。
【0014】
即ち、粘弾性材料で生じている変形状態は、加わるひずみが圧縮方向に増加する「負荷状態」と、圧縮量が徐々に減少する「復元状態」とに分けることができ、その状態によりひずみとひずみ速度の関係が異なってくる。従って、上記のように、粘弾性材料で生じているひずみ状態が負荷状態時と、除荷(あるいは復元)状態時との2つの状態を考慮し、負荷と除荷では異なったひずみ、ひずみ速度の関係を用いて粘性抵抗を決定し、ひずみがある同じ値のときでもひずみの負荷状態と除荷状態において異なる粘性抵抗値を導出することにより、実際の材料の変形状態を的確に捉えることができる。よって、材料が変形の速度、大きさ等により変化する現象を精度良く表現することができ、シミュレーションの精度をさらに向上させることができる。
【0015】
また、上記各要素毎に主ひずみの大きさ(ノルム)を計算し、シミュレーション時の前ステップのノルムと比較することによって、ノルムが増加している場合は負荷、ノルムが減少している場合は除荷と判断し、ひずみの負荷と除荷の状態を判定していることが好ましい。
例えば、一方向の主ひずみだけでは、微小な変化で負荷と除荷が様々な変化を示すが、このように、多方向の大きさにすることで、ある程度の大きさの変化でないと負荷と除荷は変化しないため、計算の安定性を高めることができる。
【0016】
本発明では、全体座標系でε、ε、ε、εxy、εyz、εzxとする要素毎に生じるひずみ、ひずみ速度の偏差成分において粘弾性特性を考慮するために、上記ひずみ、ひずみ速度を体積成分と偏差成分に分解し、主ひずみ座標系及び主ひずみ速度座標系に変換した偏差成分の主ひずみ及び主ひずみ速度を用いて主ひずみ座標系及び主ひずみ速度座標系の各軸毎にひずみの負荷状態時とひずみの除荷あるいは復元状態時とで異なった粘性抵抗を決定し、該粘性抵抗を要素毎に各軸方向3成分について可変にしている。このため、3次元方向について正確に粘弾性材料の粘弾性を考慮した解析を行うことができ、シミュレーションの精度を向上することができる。
ここで、体積成分とは、ε=(ε+ε+ε)で表される成分を指し、偏差成分とは、εxy’=εxy/2、εyz’=εyz/2、εzx’=εzx/2、ε’=ε−ε/3、ε’=ε−ε/3、ε’=ε−ε/3で表される成分を指す。なお、せん断成分を1/2倍しているのは、物理ひずみを変換しているためである。
偏差成分のみを変換しているのは、体積成分の粘弾性成分は非常に小さく、また、現状では実測で求めることができないという理由による。即ち、体積成分の粘弾性成分はほとんどないので、偏差成分のみで判断し、変換は実測の結果から1軸方向の結果しかないために、εxy、εyz、εzx成分=0になるように変換し、ε、ε、ε(主ひずみ3成分)のみの3成分について実測の結果を用いる。
なお、全体座標系とは、初期の形状、つまりモデルの形状を決定している座標系を指し、主ひずみ座標系及び主ひずみ速度座標系とは、ひずみ及びひずみ速度をせん断成分が0になるように変換した時、全体座標系を同じ角度変換した座標系を指す。
【0017】
各要素毎に生じるひずみ(ひずみ速度)は、6成分であるが、高速大変形の条件下で測定される材料物性は、1軸方向の材料物性であるために、せん断方向の材料物性の考慮が不十分である。従って、各要素毎に生じるひずみ6成分を主ひずみ座標系に変換し、主軸方向のひずみ3成分を求める。
即ち、せん断方向のひずみを0にすることで全体座標系の6成分(ε、ε、ε、εxy、εyz、εzx)を主軸方向のひずみ(主ひずみ座標系、ε、ε、ε)3成分にする。これにより、せん断方向を考慮することなく、実測結果から粘弾性材料の材料特性を評価することができる。
【0018】
具体的には、各要素毎に生じる全体座標系のひずみ及びひずみ速度を偏差成分と体積成分に分解し、偏差成分に粘弾性を考慮している。分解したひずみ、ひずみ速度の偏差6成分を、全体座標系から主ひずみ座標系及び主ひずみ速度座標系に変換し、主ひずみ座標系及び主ひずみ速度座標系での偏差主ひずみ、偏差主ひずみ速度の各3成分を用いている。これにより、主ひずみ座標系及び主ひずみ速度座標系で各軸毎に粘性抵抗を決定することができる。
【0019】
また、上記粘性抵抗から粘弾性を考慮して、偏差成分の応力とひずみを計算し、該偏差成分の応力とひずみを元の全体座標系に再変換し、体積成分の応力とひずみ及び偏差成分の応力とひずみから、要素毎の応力とひずみを算出している。このように、要素毎の応力とひずみを算出することにより、時々刻々の応力、ひずみを算出することができる。
【0020】
具体的には、下記の数式1、2より、偏差主ひずみ及び偏差主ひずみ速度を用い、各座標軸毎にひずみの負荷状態時とひずみの除荷あるいは復元状態時とで異なった粘性抵抗を決定し、該粘性抵抗から粘弾性を考慮した偏差成分の応力とひずみを算出している。
(数式1)

Figure 0003910825
ここで、E:i番目の要素の偏差主ひずみ、n:解析時間ステップ、
β:i番目の要素のβ(ここでβ=E/η、E=Eshort−Elong、η:粘性抵抗)、γ:i番目の要素の偏差主ひずみ速度を指す。
(数式2)
Figure 0003910825
ここで、σ:i番目の要素応力、F:剛性係数、Gshort:短期せん断係数、Glong:長期せん断係数を指す。
数式1、2中のβ、剛性の係数Fを実験データより算出した後、数式1により、ひずみの偏差6成分の更新を行い、数式2により応力の計算を行う。この後、主ひずみ座標系、主ひずみ速度座標系から全体座標系に変換し、下記の数式3により応力の更新を行っている。
(数式3)
Figure 0003910825
ここで、K:体積弾性係数、ε:体積ひずみ、σij:全体座標系の応力を指す。
【0021】
さらに、粘弾性材料の粘性を考慮した粘弾性モデルを用いひずみの負荷状態時とひずみの除荷あるいは復元状態時とで異なった粘性抵抗を導出し、ひずみ、ひずみ速度、ひずみの負荷状態時とひずみの除荷あるいは復元状態時とで異なる粘性抵抗の関係を製品モデルに入力しシミュレーションを行っているため、粘弾性材料の非線形性を示す物性値が、材料の変形の速度、大きさによって非線形的に変化する現象を精度良く表現することができる。かつ、シミュレーションに用いられるひずみ、ひずみ速度、応力の値は、粘弾性材料からなる製品の実使用状態を想定した条件下で測定されているため、実際の粘弾性材料の種々の変形状態に対応したシミュレーションを行うことができる。従って、材料の変形状態により、ひずみ、ひずみ速度の関係が変化し、損失係数等の材料物性が非線形性を示すような粘弾性材料からなる製品においても、シミュレーションによる性能の精度良い予測を可能にしている。
【0022】
本シミュレーション方法では、粘弾性材料からなる製品の実使用状態を想定した測定条件下で、上記粘弾性材料に生じるひずみ、ひずみ速度、応力の時々刻々の値を測定している。具体的には、実際の製品使用時に製品に外力が加わり、粘弾性材料が変形した状態を想定して測定条件を定めている。上記測定条件下で、粘弾性材料に生じるひずみ、ひずみ速度、応力の時々刻々の値を測定し、各時刻歴データを得ている。従って、上記ひずみ、ひずみ速度、応力の時刻歴データから、粘弾性材料に実使用状態を想定した外力を加わえた時の粘弾性材料の変形状態の情報を得ることができる。これにより、衝撃荷重よって高速大変形な挙動を伴う粘弾性材料についても正確に材料物性を予測することができる。
【0023】
上記ひずみ、ひずみ速度、応力の時々刻々の値は、複数の異なる測定条件下で測定するのが好ましい。製品に加わる外力の大きさを変更し、複数のパターンの測定条件を設定することにより、ひずみ、ひずみ速度、応力について種々のパターンのデータを得ることができ、シミュレーションの入力値の精度を向上することができる。また、できる限り、多くのパターンのデータを得るために、上記ひずみ、ひずみ速度、応力の値は、粘弾性材料に外力を加えて、ひずみが生じてから、ひずみがほぼゼロになるまで測定するのが好ましく、シミュレーション精度の点より測定時間の間隔は短い方が良い。
【0024】
粘弾性材料の粘性を考慮した粘弾性モデルとしては、バネとダッシュポットからなる粘弾性モデルが好ましい。このような粘弾性モデルにより、粘弾性材料の粘性を単純化することができるため、粘性が材料の変形状態に及ぼす影響を容易に考慮することができる。具体的には、マックスウェル(maxwell)モデル、フォークト(Voigt)モデルあるいは、さらに、複数のバネ、ダッシュポットの組み合わせが挙げられる。モデルの単純化という観点からは2要素モデルが好ましい。このような粘弾性材料モデルでは、ダッシュポットの粘性抵抗ηとバネの剛性(縦弾性係数Eあるいはせん断係数(横弾性係数))を可変にして用いている。せん断係数とは、縦弾性係数(ヤング率)Eとポアソン比とにより決定される物性値である。
【0025】
上記ひずみ、ひずみ速度、応力の各時刻歴データと粘弾性材料の粘性を考慮した粘弾性モデルとから、粘弾性材料のひずみの負荷状態時とひずみの除荷あるいは復元状態時とで異なった粘性抵抗の時刻歴データを導出している。具体的には、上記粘弾性モデルにより、粘弾性材料に生じるひずみ、ひずみ速度、及び粘弾性材料の縦弾性係数、ひずみの負荷状態時とひずみの除荷あるいは復元状態時とで異なった粘性抵抗の各関係を定式化している。このように定式化することで、粘性抵抗をひずみ、ひずみ速度の関数として表している。そして、上記測定により得られた粘弾性材料に生じるひずみと応力により粘弾性材料の縦弾性係数を導出し、ひずみ、ひずみ速度と共に、上記関数に代入することでひずみの負荷状態時とひずみの除荷あるいは復元状態時とで異なった粘性抵抗の値を導出している。上記測定にて、ひずみ、ひずみ速度、応力についての時刻歴データを得ているため、粘性抵抗についても同様に、対応する時刻歴データを得ることができる。このように、粘性抵抗は、ひずみ、ひずみ速度の値によって決定され、その時間変化に対応して粘性抵抗も変化する。
【0026】
上記粘弾性材料からなる製品を解析対象の製品モデルとして設定し、該製品モデルを含み、速度、拘束条件等を含んだ計算用入力データ(あるいはinput data)に上記ひずみ、ひずみ速度、ひずみの負荷状態時とひずみの除荷あるいは復元状態時とで異なった粘性抵抗の関係を入力し、上記ひずみ、ひずみ速度の違いによる粘性抵抗の変化を考慮したシミュレーションを行っている。
シミュレーション計算時に製品モデルに上記ひずみ、ひずみ速度、ひずみの負荷状態時とひずみの除荷あるいは復元状態時とで異なった粘性抵抗の関係を入力しているが、具体的には、ひずみとひずみ速度の関係、及び、ひずみと粘性抵抗の関係の各2次元データを入力し、演算することができる。また、ひずみとひずみ速度と粘性抵抗の関係を3次元データとし、粘性抵抗を関数化したデータとして入力し、演算することもできる。
【0027】
ひずみとひずみ速度の関係、及び、ひずみと粘性抵抗の関係の各2次元データは、上記各関係を用いて、ひずみ、ひずみ速度とそれに対応するひずみの負荷状態時とひずみの除荷あるいは復元状態時とで異なった粘性抵抗を入力データとして書き込むことになる。具体的には、複数の測定条件にて、ひずみ、ひずみ速度等の測定を行い、測定条件の異なる各パターンについて、ひずみ、ひずみ速度の時系列データから、各々ひずみとひずみ速度の対応関係を記録する。各々の曲線に対応する粘性抵抗の値も記録する。これらひずみ、ひずみ速度、粘性抵抗の3つの対応を整理することで、あるひずみ、ひずみ速度の時の粘性抵抗を精度良く導出し、計算を行っている。
【0028】
異なるひずみ、ひずみ速度の条件下での測定数が多いほど精度良く粘弾性材料の物性を実現できるため、上記したように複数の異なった測定条件にてひずみ、ひずみ速度等の測定を行うのが好ましい。ただし、測定データが多いほど、シミュレーションの際、計算時間を要するため、定められた測定条件で測定した測定データと同一でないひずみ、ひずみ速度の場合には、補間を用いて粘性抵抗を算出するのが好ましい。上記補間は、例えば、近接するひずみ、ひずみ速度の条件下での粘性抵抗2値からの1次補間、あるいは測定した全ての条件下での個々の測定値を用いての補間等の種々の方法により行うことができる。このような補間作業を行うことにより、測定条件の違いによる材料に生じるひずみとひずみ速度の変化に対応して粘性抵抗を算出することができる。
【0029】
上記粘弾性材料は、材料物性が非線形性を示す材料であっても本シミュレーション方法によれば、実使用状態を想定した材料物性及び変形挙動を正確にシミュレーションすることができる。粘弾性材料の材料物性が、一次直線的でなく、材料の変形速度や変形の大きさによって変化を伴う強い非線形性を示す材料においても、上記のように、粘性抵抗をひずみ、ひずみ速度の2つの値によって決定するような粘弾性モデルを用い粘性抵抗を算出することで、非線形性を考慮することができる。本シミュレーション方法では、特に、損失係数が強い非線形性を示す材料において、精度良くシミュレーションを行い、実使用状態を想定した製品の性能を予測することができる。
【0030】
本発明のシミュレーション方法は、以下に示すように有限要素法により行っている。
上記シミュレーションを有限要素法解析により行う際は、上記製品モデルに多数の節点と要素を設定している。即ち、有限要素法で粘弾性材料で構成される製品をシミュレーションして材料物性を予測する際に、粘弾性モデルのバネの剛性(縦弾性係数あるいはせん断係数)を決定し、粘性を示す粘性抵抗を上記に示すように各要素ごとにその要素で生じているひずみ、ひずみ速度によって決定することで各要素ごとに適正なひずみ、ひずみ速度の条件下の材料物性を表すことができる。なお、縦弾性係数を入力する代わりに、ポアソン比との関係から、せん断係数で入力しても構わない。縦弾性係数で入力するか、せん断係数で入力するかは、有限要素法プログラムの仕様により変更できる。
【0031】
上記ひずみ、ひずみ速度、応力をスプリットホプキンソン棒試験機により測定している。スプリットホプキンソン棒試験機を用いると、試験片に高速でかつ大変形のひずみを与えることができ、1万分あるいは1千分の数秒という高速で且つ変形量も数十%という、高速大変形条件下での粘弾性材料のひずみ、ひずみ速度、応力の時系列データを得ることができる。スプリットホプキンソン棒試験機の測定条件を、製品に衝撃荷重が加わった時の粘弾性材料に生じるひずみ、ひずみ速度と同等の条件とすれば、外力により生じた粘弾性材料の高速大変形状態等、種々の状態に対応した粘弾性材料の材料物性を得ることができる。従って、このスプリットホプキンソン棒試験機で測定された材料物性を用いることで、シミュレーションの精度を向上することができる。
【0032】
また、スプリットホプキンソン棒試験機は、試験片に衝撃を加える打撃棒の衝突速度を変更するだけで、様々なひずみ、ひずみ速度の領域での試験片の材料物性を測定することができるため、容易に、種々のひずみ、ひずみ速度のパターンにおいて、材料物性を得ることができる。
【0033】
スプリットホプキンソン棒試験機は、元来、金属材料の衝撃挙動の評価用に用いられていたが、本発明では、スプリットホプキンソン棒試験機を、粘性を持つ粘弾性材料の評価用に改良して用いている。スプリットホプキンソン棒試験機の測定方法等については後述する。
【0034】
なお、試験片に高速でかつ大変形のひずみを与えることができ、製品が実際に使用される条件下でのひずみ、ひずみ速度の条件下で、材料物性を測定することができれば、上記スプリットホプキンソン棒以外の測定方法により材料物性を求めても良いことは言うまでもない。
【0035】
上記粘弾性材料からなる製品の実使用状態を想定した測定条件下でのひずみ、ひずみ速度、応力の測定時に、粘弾性材料に生じる、上記ひずみの最大値は0.05〜0.50の範囲であり、あるいは/及び、上記ひずみ速度の最大値は500/s〜10000/s、好ましくは500/s〜5000/sの範囲であるのが良い。上記ひずみの最大値と、ひずみ速度の最大値の範囲は、粘弾性材料の高速大変形時に生じるひずみ、ひずみ速度の条件であるため、高速大変形時の製品の性能を予測するには、この条件下のひずみ、ひずみ速度、応力の3つの時系列データを用いることが好ましい。
【0036】
上記粘弾性材料はゴルフボール用材料とし、上記製品モデルをゴルフボールとしている。ゴルフボールは、粘弾性材料から構成され、実使用時に衝撃荷重等の外力を受け高速変形、あるいは大変形を強いられる製品であり、その高速、大変形状態がゴルフボール自体の性能に大きく影響を及ぼす。このため、本発明のシミュレーション方法による解析は、ゴルフボールの性能予測に非常に有用であり、試作をせずに、ゴルフボールの性能を精度良く予測することができるため、ゴルフボールの設計の効率化を図ることができる。
【0037】
上記ゴルフボールと、ゴルフクラブヘッドを想定した打撃物との衝突現象をシミュレーションし、該衝突時のゴルフボールの挙動を予測している。実際にゴルフクラブヘッドにより、ゴルフボールが打撃されたときと同等のひずみ、ひずみ速度、応力の条件下で、ゴルフボールを構成する粘弾性材料の材料物性を予測することができるため、ゴルフボールの反発係数やゴルフボール打撃時のゴルフボールの変形状態等、ゴルフボールの挙動を予測することができる。
【0038】
上記ゴルフボールは、架橋ゴム層等の単体からなる所謂1ピースゴルフボールであってもよく、架橋ゴム層等のコアにカバーが被覆された所謂2ピースゴルフボールでもよく、また、3層以上から構成されていてる所謂マルチピースゴルフボールであってもよく、粘弾性材料からなるあらゆる構造のゴルフボールに適応可能である。
【0039】
上記粘弾性材料としては、粘弾性を有するあらゆる材料が含まれる。例えば、熱可塑性樹脂、熱硬化性樹脂、各種エラストマー、ゴム等が挙げられ、これらの単体、あるいは混合物を用いることができる。また、これらの単体、混合物には、着色剤、劣化防止剤、架橋剤等の各種添加剤が、必要に応じ配合されていてもよい。製品の実使用条件下で、ひずみ、ひずみ速度、応力を測定できる材料であれば、あらゆる材料に適応可能である。
【0040】
具体的には、ゴルフボールの構成材料として用いられるアイオノマー樹脂等の合成樹脂、ポリブタジエン(ブタジエンゴム)、天然ゴム、ポリイソプレン、スチレン−ブタジエン共重合体、エチレン−プロピレン−ジエン共重合体(EPDM)、ウレタンゴム等が挙げられる。
【0041】
粘弾性材料からなる製品としては、ゴルフボール以外に、プリンタ等の印刷機に用いられるゴムローラ、ゴムベルト、あるいはタイヤ、その他粘弾性材料からなるテニス用、ゴルフ用等のスポーツ用品等が挙げられる。粘弾性材料からなる製品としては、製品の少なくとも一部に粘弾性材料が使用されていれば良く、金属材料等の他の材料との複合成形品等であっても良く、製品の粘弾性材料に該当する部分の性能を予測することもできる。特に、使用時において、衝撃荷重を受け、高速で大きな変形量を示すような製品に好適であり、本シミュレーション方法により製品の性能や動的挙動を精度良く予測することができる。
【0042】
また、本シミュレーション方法においては、上記製品モデルを含む計算用入力データに上記ひずみ、ひずみ速度、縦弾性係数あるいはせん断係数を入力し、上記ひずみ、ひずみ速度の違いによる縦弾性係数あるいはせん断係数の変化を考慮することにより、縦弾性係数が材料の変形状態(ひずみ、ひずみ速度の大きさ等)により変化する粘弾性材料についても正確な解析を行うことができ、シミュレーションの精度をさらに向上させることができる。
また、場合によっては、応力の負荷状態と除荷状態で場合分けすることもできる。
【0043】
【発明の実施の形態】
以下、本発明の実施形態を図面を参照して説明する。
本発明の第一実施形態のシミュレーション方法では、非線形性を示す粘弾性材料として、ゴルフボールの構成材料であるブタジエンゴムを主成分とする材料を用いている。上記ブタジエンゴムを主成分とする材料を試験片として、粘弾性材料測定用に改良したスプリットホプキンソン棒試験機によって、ゴルフボールの実使用状態を想定した高速で大きな変形時の測定条件下で、上記ブタジエンゴムを主成分とする材料に生じるひずみ、ひずみ速度、応力の時々刻々の値を測定している。なお、スプリットホプキンソン棒試験機による測定方法については、後述する。
【0044】
スプリットホプキンソン棒試験機は打撃棒の衝突速度を変更することで様々なひずみ、ひずみ速度の領域での材料物性を測定することができる。本実施形態では、4パターンの衝突速度(7m/s、14m/s、20m/s、25m/s)を採用し、4つの測定条件下で測定を行い、それぞれの衝突パターンにてひずみ、ひずみ速度、応力の時刻歴データを得ている。上記ひずみε、ひずみ速度ε’、応力σの各時刻歴データを、それぞれ図1、図2、図3に示す。
【0045】
スプリットホプキンソン棒試験機により打撃棒を衝突させた後、図1に示すように、ブタジエンゴムを主成分とする材料にはひずみが生じており、一定時間まではひずみの値は増加し、その後、なだらかに減少している。また、図2に示すように、ひずみ速度はひずみが最大値に達するまでは正の値を示し、ひずみの最大値以降は、負の値を示している。図3に示すように、応力の値も同様に、時間と共に変化している。
【0046】
上記ひずみの時刻歴データと応力の時刻歴データからひずみ−応力線図を描くと図4のようになる。図4のひずみ−応力線図において、ひずみの値が最大時のひずみεmaxと、その時の応力σaの値を用い、下記の数式4により、各衝突パターン毎に試験片の縦弾性係数Eを算出している。
(数式4)
E=σ/εmax −−−(4)
なお、本実施形態で用いたブタジエンゴムを主成分とする材料では、スプリットホプキンソン棒試験機において打撃棒の衝突速度を変更しても、各衝突パターンでの縦弾性係数はほぼ同じ値を示している。
【0047】
粘弾性材料からなる製品について、粘性を考慮したシミュレーションを行うために、粘弾性材料の粘性を考慮した粘弾性モデルを設定している。具体的に、本実施形態では、図5に示すような、バネとダッシュポットからなる粘弾性モデルとして、基本的な2要素のVoigtモデルを用いている。即ち、上記粘弾性材料モデルでは、ダッシュポットの粘性抵抗ηとバネの剛性(縦弾性係数Eあるいはせん断係数)を可変にして用いている。
【0048】
図5に示すように2要素のフォークト(Voigt)モデルとした粘弾性モデルにて、バネで生じる応力をσ, ダッシュポットで生じる応力をσとすると全体で生じる応力σは以下のような数式5で表すことができる。
(数式5)
σ=σ+σ
=Eε+ηε’ −−−(5)
【0049】
従って、上記数式5より、このような粘弾性モデルにて粘性抵抗ηは、以下の数式6で表すことができる。
(数式6)
η=(σ−Eε)/ε’ −−−(6)
【0050】
従って、図1、図2、図3に示すひずみε、ひずみ速度ε’、応力σの時刻歴データと上記の数式3により、ひずみε、ひずみ速度ε’に対応する粘性抵抗ηを時々刻々算出できる。ここで、上記したように、上記ひずみ、ひずみ速度、応力の値は、時間の経過と共に変化しており、粘弾性材料で生じているひずみ状態は、加わるひずみが圧縮方向に増加する「負荷状態」と、圧縮量が徐々に減少する「復元状態」とに分けることができる。よって、粘性抵抗を算出する際にも、粘弾性材料で生じているひずみ状態が負荷状態時と、除荷(あるいは復元)状態時との2つに分けて算出し、図6に示すように、各衝突速度において粘性抵抗の時刻歴データを得ている。
【0051】
このように、上記高速大変形条件下での測定により得られたひずみと応力より、粘弾性材料の縦弾性係数Eを決定した後、粘性を考慮した粘弾性モデルから導出した上記数式6により、負荷状態時と、除荷(あるいは復元)状態時とに分けて粘性抵抗ηを算出し、このひずみの負荷状態時とひずみの除荷あるいは復元状態時とで異なった粘性抵抗ηをシミュレーションに適用する。
【0052】
上記ひずみの時刻歴データと粘性抵抗の時刻歴データより、ひずみと粘性抵抗との関係を図7に示す。図7に示すように、粘性抵抗の値は、ひずみの負荷状態時と除荷(あるいは復元)状態時とでは、ひずみが同じ値のときでも異なった値を示している。また、ひずみの時刻歴データとひずみ速度の時刻歴データより、ひずみとひずみ速度との関係を図8に示す。このように、粘性抵抗は、ひずみ、ひずみ速度によって変化し、かつ、ひずみの負荷状態時と除荷(あるいは復元)状態時とで異なる値を示すため、粘性抵抗はひずみ、ひずみ速度によって決定することができ、粘性抵抗をひずみの負荷状態時と除荷状態時とで場合分けして、ひずみ、ひずみ速度の関数としている。なお、図8において、4つの衝突速度の各データは、ひずみ速度の負の領域では、表示上一直線に重なっているように見えるが、実際の値は、各衝突速度毎にそれぞれ異なるひずみ、ひずみ速度の値となっている。
【0053】
図7に示すひずみと粘性抵抗の関係と、図8に示すひずみとひずみ速度の関係とを、後述するように各衝突パターンにおいて、それぞれ、入力データとして書き込み、有限要素法解析により上記ひずみ、ひずみ速度の違いによる粘性抵抗の変化、即ち、ひずみの負荷状態時と除荷(あるいは復元)状態時との違いによる粘性抵抗の変化を考慮したシミュレーションを行っている。
【0054】
本実施形態では、粘弾性材料からなる製品モデルとしてゴルフボールモデルを設定し、ゴルフクラブヘッド(打撃物)のゴルフボールへの衝突を想定したシミュレーションを行っている。図9に示すように、ゴルフボールモデル10は、ブタジエンゴムを主成分とする1ピースボールを想定し、直径42.8mmの球状としている。
【0055】
有限要素法解析を行うにあたり、製品モデルにおいて初期条件を設定している。即ち、ゴルフボールモデル10の大きさ、形状、構造等の初期条件を設定すると共に、ゴルフボールモデル10を多数の要素11にメッシュ分割し、多数の節点12を得ている。ゴルフボールモデルの全ての範囲の要素の数は、全体モデルで1000〜100000が好ましく、さらに好ましくは、2500〜20000が良い。なお、上限の値は、現段階での計算機の能力を鑑みてのものであり、今後計算機の能力が向上するにつれ、解析の時間が短縮されるため上限の範囲は、変わることが容易に想像できる。
【0056】
上記設定条件に基づいて、シミュレーション計算時に、上記ひずみ、ひずみ速度、粘性抵抗の関係とゴルフボールモデル10のデータを入力データに書き込む。計算が実行されると時々刻々各要素に関して妥当な粘性抵抗を算出し、その粘性抵抗を用い式の計算を行っている。本シミュレーションでは、ひずみの負荷状態時と除荷(あるいは復元)状態時とで場合分けされたひずみとひずみ速度の関係、及び、ひずみと粘性抵抗の関係の各2次元データを入力し演算している。測定条件の異なる4つの各パターンについてのひずみとひずみ速度の時系列データから、各々ひずみとひずみ速度の対応関係を記録し、ひずみ、ひずみ速度と、それに対応する粘性抵抗を入力データとして書き込むことになる。4つの測定条件により測定されたひずみ、ひずみ速度の測定データと同一でないひずみ、ひずみ速度の場合には、近接するひずみ、ひずみ速度の条件下での粘性抵抗2値を用いて1次補間を行っている。
【0057】
さらに具体的には、一つの要素に着目したとき、その要素に生じているひずみの大きさ(ノルム)を計算し、シミュレーション時の前ステップのノルムと比較することによって、ノルムが増加している場合は負荷、ノルムが減少している場合は除荷と判断し、ひずみの負荷と除荷の状態を判定し、材料状態を考慮したひずみとひずみ速度の情報を得ている。次に、測定から得られた各衝突ケース(測定ケース)のひずみ、ひずみ速度の関係から、各ケースでの同値のひずみが生じているときのひずみ速度の値を参照して、着目する要素のひずみ速度の値に近いひずみ速度の値を持つ2つのケースを探索する。この2つのケースに対応するひずみが同値のときの,ひずみ速度と粘性抵抗の値を用いて、補間により、その要素で生じるひずみ、ひずみ速度に対応した適切な粘性抵抗の値を算出することができる。本実施形態では、このように単純な一次補間による補間を行っている。このように、各要素ごとにその要素で生じているひずみ、ひずみ速度によって負荷と除荷では異なった粘性抵抗を決定することで、各要素ごとに適正なひずみ、ひずみ速度の条件下の材料物性を表している。
【0058】
また、上記要素毎に生じる全体座標系のひずみ及びひずみ速度を偏差成分と体積成分に分解し、上記偏差成分のひずみ及びひずみ速度を、全体座標系から主ひずみ座標系及び主ひずみ速度座標系に変換し、変換した偏差主ひずみ及び偏差主ひずみ速度を用い、上述したように図7、8より各座標軸毎に、ひずみが同じ値でもひずみの負荷状態時とひずみの除荷あるいは復元状態時とで異なった粘性抵抗を決定し、該数式1、2、3により粘性抵抗から粘弾性を考慮した偏差成分の応力とひずみを算出している。その後、上記偏差成分の応力とひずみを全体座標系に再変換し、上記体積成分の応力とひずみ、及び偏差成分の応力とひずみから、要素毎の応力とひずみを算出している。
【0059】
即ち、各要素毎に生じるひずみ6成分を主ひずみ座標系に変換し、主軸方向のひずみ3成分を求めている。各要素毎に生じる全体座標系のひずみ及びひずみ速度を偏差成分と体積成分に分解し、偏差成分に粘弾性を考慮している。分解したひずみ、ひずみ速度の偏差6成分を、全体座標系から主ひずみ座標系及び主ひずみ速度座標系に変換し、主ひずみ座標系及び主ひずみ速度座標系での偏差主ひずみ、偏差主ひずみ速度の各3成分を用い、主ひずみ座標系及び主ひずみ速度座標系で各軸毎に、ひずみが同じ値でもひずみの負荷状態時とひずみの除荷あるいは復元状態時とで異なった粘性抵抗を決定している。
【0060】
本実施形態では、具体的に、図10(A)(B)(C)に示すように、ゴルフボールモデル10に、打撃物として200g(ゴルフクラブヘッドの重さと同一)の円筒形状のアルミニウム製中空棒モデル20を速度45m/sで衝突させた時のゴルフボールモデル10の材料物性を有限要素法により解析し、シミュレーションを行っている。これにより、打撃物の衝突時から所定時間後までの所定時間におけるゴルフボールモデル10の各要素11に発生するひずみ量を演算している。
【0061】
1要素に関して構成する各節点に質量を分配し、各節点を質点と置き換え、各節点の速度を質点の持つ速度として合計を割って、全体の速度としている。即ち、インパクト後のボールの速度Vbi(i=x、y、z)は、数式7により以下のように算出する。ボールの全運動量は全質点の運動量の和と考え、全運動量を総重量で割ったものをボールの飛び出し速度Vbiと定義する。
【0062】
(数式7)
Figure 0003910825
ここで、N、Mは、ボールの節点数、総質量、Vni、Mnはn番目の並進速度、節点を含む要素の質量をその要素に含まれる節点の数で割ったものである。
【0063】
なお、アルミニウム製中空棒モデル20の円形面20aを衝突面とし、衝突面はフラットであり、ゴルフボールモデル10とフラットな正面衝突としている。また、アルミニウム製中空棒モデル20の円形面20aの中心点20bが、最初のゴルフボールモデル10との衝突点となるように衝突させている。
【0064】
上記した方法により、衝突前後の上記アルミニウム製中空棒モデル20及びゴルフボール10の速度を算出し、それぞれの速度及び重量からゴルフボール10の反発係数を算出し予測している。
【0065】
このように、上記ひずみ、ひずみ速度を体積成分と偏差成分に分解し、主ひずみ座標系及び主ひずみ速度座標系に変換した偏差成分の主ひずみ及び主ひずみ速度を用いて主ひずみ座標系及び主ひずみ速度座標系の各軸毎に、ひずみが同じ値でもひずみの負荷状態時とひずみの除荷あるいは復元状態時とで異なった粘性抵抗を決定し、上記粘性抵抗を要素毎に各軸方向3成分について可変にしている。このため、3次元方向について正確に粘弾性材料の粘弾性を考慮した解析を行うことができ、シミュレーションの精度を向上することができる。
【0066】
また、ひずみの負荷状態時と除荷(あるいは復元)状態時とで場合分けされた粘性抵抗をひずみ、ひずみ速度の関数として決定し、上記ひずみ、ひずみ速度、粘性抵抗の3者の関係をゴルフボールモデルに入力して、有限要素法解析によりシミュレーションしている。このため、各要素の時々刻々のひずみ、ひずみ速度からそれに相当する粘性抵抗を算出でき、該粘性抵抗値はひずみの負荷状態時と除荷(あるいは復元)状態時との違いが考慮されており、高速大変形条件下においても強い非線形性を持つ粘弾性材料の特性を非常に精度良く予測することができる。
【0067】
従って、実際のゴルフクラブヘッドによる打撃時と同等の条件下における、粘弾性材料からなるゴルフボールモデルの反発係数等の性能、動的挙動等を、試作をすることなく正確に把握することができる。
【0068】
また、図11示すように、応力−ひずみ曲線から、下記の数式8を用いて、位相角δが算出され、この位相角δより、損失係数(tanδ)を算出することもできる。
(数式8)
δ=sin−1((σ−σ)/σmax) −−−(8)
【0069】
(スプリットホプキンソン棒試験機による材料物性の測定)
図12は、粘弾性材料測定用に改良したスプリットホプキンソン棒試験機が示された模式的正面図である。
【0070】
図12に示すスプリットホプキンソン棒試験機は、打撃棒51、入力棒53及び出力棒55を備えており、これらは直線上に配置されている。入力棒53には、第一ひずみゲージ57及び第二ひずみゲージ59が取り付けられている。出力棒55には、第三ひずみゲージ61及び第四ひずみゲージ63が取り付けられている。入力棒53の後端53aと出力棒55の前端55aとの間には、円柱状の試験片70が挟持されている。
【0071】
試験片70は、測定すべき粘弾性材料を試験片の形状に成形したものでもよく、また、粘弾性材料を成形した製品状態から切り出されたものであってもよい。本実施形態では、試験片70の長さ(すなわち入力棒53の後端53aと出力棒55の前端55aとの距離)は4mm、試験片70の断面直径は18mmとしている。
【0072】
打撃棒51、入力棒53及び出力棒55はポリメチルメタアクリレート製の円柱であり、断面直径は20mm、縦弾性係数は5300MPa、比重は1.19である。打撃棒51の長さは、100mmである。入力棒53及び出力棒55(以下、この入力棒53と出力棒55とは、「応力棒」とも称される)の長さは、2000mmである。
第一ひずみゲージ57は入力棒53の後端53aから900mmの位置に取り付けられており、第二ひずみゲージ59は入力棒53の後端53aから600mmの位置に取り付けられている。また、第三ひずみゲージ61は出力棒55の前端55aから300mmの位置に取り付けられており、第四ひずみゲージ63は出力棒55の前端55aから600mmの位置に取り付けられている。
【0073】
このように、図12のスプリットホプキンソン棒試験機は、打撃棒51、入力棒53及び出力棒55がポリメチルメタアクリレートからなる合成樹脂製であり、打撃棒51及び入力棒53が2000mmと大きく、しかも第一ひずみゲージ57と入力棒53の後端53aとの距離及び第二ひずみゲージ59と入力棒53の後端53aとの距離が大きいので、ゴルフボールに用いられるような架橋ゴム等の粘弾性を有する材料のひずみ、ひずみ速度、応力の測定に適している。
【0074】
上記第一ひずみゲージ57、第二ひずみゲージ59、第三ひずみゲージ61、第四ひずみゲージ63として単軸プラスチック用ひずみゲージを用い、本実施形態では、株式会社共和電業製のKFP−5−350−C1−65を用い、入力棒53、出力棒55の上記した位置に貼着している。これら第一ひずみゲージ57乃至第4ひずみゲージ63の入力棒53及び出力棒57への取付位置は、長さ方向において同一線上としている。
【0075】
このスプリットホプキンソン棒試験機によるひずみ、ひずみ速度、応力の測定では、まず、打撃棒51が入力棒53の前端53aに、所定の速度で衝突する。これによって入射ひずみ波が生じ、この入射ひずみ波は入力棒53の後端53aに向かって進む。この入射ひずみ波の一部は、入力棒53の後端53aにおいて反射し、反射ひずみ波となって入力棒53の前端53bに向かって進む。入射ひずみ波の一部は、入力棒53の後端53aから試験片70を透過し、さらに出力棒55に伝播して透過ひずみ波となり、出力棒55の後端55bに向かって進む。
【0076】
入射ひずみ波は、第一ひずみゲージ57及び第二ひずみゲージ59によって実測される。実測された入射ひずみ波は、ローパスフィルターに通され、10kHz以上の高周波が除去される。さらに、入射ひずみ波の時刻歴は、そのベースライン値をゼロとするゼロ補正が施される。こうして得られた第一ひずみゲージ57及び第二ひずみゲージ59における時間軸ひずみのそれぞれがフーリエ変換され、周波数軸ひずみが求められる。この第一ひずみゲージ57及び第二ひずみゲージ59における周波数軸ひずみから、伝達関数が導出される。第一ひずみゲージ57と入力棒53の後端53aとの距離X1と、第二ひずみゲージ59と入力棒53の後端53aとの距離X2との比(X1:X2)が考慮されつつ、上記伝達関数に基づいて、入力棒53の後端53aにおける周波数軸ひずみが推定される。この周波数軸ひずみがフーリエ逆変換されることにより、入力棒53の後端53aにおける入射ひずみ波の時間軸ひずみ(ひずみの時刻歴)εiが得られる。
【0077】
同様に、入力棒53の後端53aで反射して前端53bに向かう反射ひずみ波が第二ひずみゲージ59及び第一ひずみゲージ57によって実測される。この実測された反射ひずみ波から、入力棒53の後端53aにおける反射ひずみ波の時間軸ひずみ(ひずみの時刻歴)εrが得られる。
【0078】
また、出力棒55の第三ひずみゲージ61及び第四ひずみゲージ63によって、試験片70を経て出力棒55に伝播される透過ひずみ波を実測し、この実測した透過ひずみ波から、出力棒55の前端55aにおける透過ひずみ波の時間軸ひずみ(ひずみの時刻歴)εtが得られる。
【0079】
こうして得られたεi、εr及びεtから、下記数式9によって、試験片70のひずみ速度ε’が算出される。
(数式9)
ε’=(C/L)・(εi−εr−εt)
=((E/ρ)1/2/L)・(εi−εr−εt) −−−(9)
(数式9において、Cは入力棒および出力棒中(応力棒)のひずみ波の伝播速度(m/s)を表し、Lは試験片の長さ(m)を表し、Eは応力棒の縦弾性係数(N/m)を表し、ρは応力棒の密度(kg/m)を表す)
【0080】
また、εi、εr及びεtから、下記数式10によって試験片70のひずみεが算出される。
【0081】
(数式10)
Figure 0003910825
(数式10において、Cは入力棒および出力棒中(応力棒)のひずみ波の伝播速度(m/s)を表し、Lは試験片の長さ(m)を表し、Eは応力棒の縦弾性係数(N/m)を表し、ρは応力棒の密度(kg/m)を表す)
【0082】
さらに、εi、εr及びεtから、下記数式11によって試験片70の応力σが算出される。
(数式11)
σ=(E・A/(2As))・(εi+εr+εt)
=(E・D/(2(Ds)))・(εi+εr+εt)−−−(11)
(数式11において、Eは入力棒および出力棒からなる応力棒の縦弾性係数(N/m)を表し、Aは上記応力棒の断面積(m)を表し、Asは試験片の断面積(m)を表し、Dは応力棒の直径(m)を表し、Dsは試験片の直径(m)を表す)
【0083】
こうして得られた試験片70のひずみ時刻歴を、図13のグラフに示す。図13に示すように、曲線は、ピークP以降しばらくはなだらかであるが、その後、凹凸状となる。ピークP以降のなだらかな段階での点Sを選択し、この点Sにおける曲線に対する接線を画き、この接線と時間軸との交点から緩和時間λを導出し、下記数式12によって求められる曲線を点S以降の曲線とすることによって、ひずみ時刻歴全体をなだらかな曲線(図13中に点線で示す)とすることができる。これにより、最終的に得られる材料物性へのノイズの影響を除去することができる。
(数式12)
ε(t)=ε・e−t/λ −−−(12)
(数式12において、εは接点におけるひずみを表す)
【0084】
同様に、下記数式13によって、応力時刻歴全体をなだらかな曲線とすることができ、これによって最終的に得られる材料物性へのノイズの影響を除去することができる。
(数式13)
σ(t)=σ・e−t/λ −−−(13)
(数式(13)において、σは接点における応力を表す)
【0085】
上記のようにして、試験片70のひずみ時刻歴及び応力時刻歴については補正が行われている。
【0086】
以上の方法により、スプリットホプキンソン棒試験機にて、高速大変形時のひずみの時刻歴データ、ひずみ速度の時刻歴データ、応力の時刻歴データを得ている。
【0087】
以下、本発明の実施例及び比較例について詳述する。
まず、ウレタンゴムを主成分とする材料を用いてゴルフボールを作製した。ウレタンを主成分とする材料を160℃で30分、直径42.8mmの金型で圧縮成型し、1ピースボールとした。
【0088】
上記ウレタンゴムを主成分とする材料のひずみ、ひずみ速度、応力を上記したスプリットホプキンソン棒試験機により、打撃棒の衝突速度が7m/s、14m/s、20m/s、25m/sの4パターンについて測定した(室温23℃、相対湿度50%)。
測定時の最大ひずみと最大ひずみ速度を以下に示す。
衝突速度7m/s(最大ひずみ0.12、最大ひずみ速度1378/s)
衝突速度14m/s(最大ひずみ0.24、最大ひずみ速度2703/s)
衝突速度20m/s(最大ひずみ0.35、最大ひずみ速度3898/s)
衝突速度25m/s(最大ひずみ0.43、最大ひずみ速度4716/s)
4パターンの衝突速度で測定しているため、各衝突速度における位相角δの値を下記の表1に示す。
【0089】
【表1】
Figure 0003910825
【0090】
(実施例)
ウレタンゴムを主成分とする材料を試験片としてスプリットホプキンソン棒試験機により測定した、ひずみ、ひずみ速度、応力の各時刻歴データと、第1実施形態と同様の粘弾性モデルとを用い、粘性抵抗を考慮してシミュレーションを行った。上記第1実施形態のシミュレーション方法において、さらに、上記製品モデルに上記ひずみ、ひずみ速度、縦弾性係数の関係を入力することで、上記ひずみ、ひずみ速度の違いによる縦弾性係数の変化を考慮した。使用ソフトは、日本総合研究所(株)製、汎用ソフトLs−dyna940とした。
実施例のシミュレーションにより予測された位相角δを上記表1に示す。
【0091】
(比較例)
上記実施例と同様のウレタンゴムを主成分とする材料を試験片としてスプリットホプキンソン棒試験機により測定した、ひずみ、ひずみ速度、応力の各時刻歴データを用いた。縦弾性係数は可変とし、負荷状態時のひずみ、ひずみ速度と負荷の粘性抵抗を利用したシミュレーションを行った。上記した点以外は実施例と同様とした。
比較例のシミュレーションにより予測された位相角δを上記表1に示す。
【0092】
有限要素法解析により、上記材料ウレタンゴムを主成分とする材料よりなるゴルフボールに、200g(ヘッドの重さと同一)のアルミニウム製の中空棒モデルを初速度35、40、45m/sの3通りの速度で衝突させた時のゴルフボールの性能、材料の変形状態をシミュレートし、上記した方法でゴルフボールの解析での反発係数を算出した。ここで得られた上記ウレタンゴムを主成分とする材料よりなるゴルフボールの実施例及び比較例の解析での反発係数を、それぞれ下記の表2に示す。
【0093】
【表2】
Figure 0003910825
【0094】
また、上記ウレタンゴムを主成分とする材料で成型されたゴルフボールの実物を使用した実験により、下記の方法でゴルフボールの実物での反発係数の値を測定した。
上記実施例及び比較例の解析での反発係数と下記の実物を使用した実験での反発係数との差(%)を上記した表2に示す。
【0095】
(ゴルフボールの実物を用いた実験による反発係数の測定)
ゴルフボールの反発係数を測定する方法として、室温23℃の条件下、上記材料よりなるゴルフボールに、ゴルフクラブのヘッドの代用として、200g(ヘッドの重さと同一)のアルミニウム製の中空棒を速度35、40、45m/sの3通りの速度で衝突させた。衝突前後の上記中空棒及びゴルフボールの速度を測定し、それぞれの速度及び重量からゴルフボールの反発係数を算出した。
なお、アルミニウム製の中空棒の衝突面は、フラットであり、ゴルフボールとフラットな正面衝突とした。ゴルフクラブヘッドのように衝突面に角度がないため、衝突時にゴルフボールが回転しないので、ゴルフボールの反発性能のみを評価可能とした。
【0096】
表1に示すように、実施例の解析での位相角δは、各衝突速度において、実験結果の位相角δとほぼ同一の値となった。一方、比較例の解析での位相角δは、負荷と負荷の曲線を用いてシミュレーションしているために、応力−ひずみ曲線においてエネルギーロスが実施例に比べて小さくなり、各衝突速度において比較例の位相角δの値は実施例より小さくなり、実験結果に比べその値も小さく、その差も大きかった。従って、実施例は精度良く実際の実験結果をシミュレーションできていることが確認できた。
【0097】
また、表2に示すように、実施例の解析での反発係数の値は、実際にゴルフボールを試作して、実物により求めた実験での反発係数の値と、いずれの中空棒速度の場合においても、よく一致していた。実施例と実物実験結果との差は−4.82%〜−5.98%であり、実物の実験結果を正確にシミュレートしていることが確認できた。一方、比較例の解析での反発係数の値は、位相角δの値が小さいために、反発係数は実験値に対して大きくなり、実施例に比べてその差も大きかった。具体的には、比較例と実物実験結果との差は16.1%〜22.2%であり、実験結果と大きく異なっていた。
【0098】
以上より、ひずみの負荷状態時とひずみの除荷あるいは復元状態時とで場合分けして、ひずみが同じ値でもひずみの負荷状態時とひずみの除荷あるいは復元状態時とでの粘性抵抗の違いを考慮した本発明のシミュレーションを行うことで、材料物性が非線形性を示す粘弾性材料からなる製品の実使用条件下での性能を、精度良くシミュレーションにより予測できることが確認できた。
【0099】
【発明の効果】
以上の説明より明らかなように、本発明によれば、要素毎に生じるひずみ、ひずみ速度の偏差成分において粘弾性特性を考慮するために、上記ひずみ、ひずみ速度を体積成分と偏差成分に分解し、主ひずみ座標系及び主ひずみ速度座標系に変換した偏差成分の主ひずみ及び主ひずみ速度を用いて主ひずみ座標系及び主ひずみ速度座標系の各軸毎に粘性抵抗を決定し、粘性抵抗を要素毎に各軸方向3成分について可変にしている。このため、3次元方向について正確に粘弾性材料の粘弾性を考慮した解析を行うことができ、シミュレーションの精度を向上することができる。
【0100】
また、粘弾性材料の粘性を考慮した粘弾性モデルを用い、ひずみの負荷状態時とひずみの除荷あるいは復元状態時とで場合分けして粘性抵抗を導出し、ひずみ、ひずみ速度、ひずみの負荷状態時とひずみの除荷あるいは復元状態時とで異なった粘性抵抗の関係を製品モデルに入力しシミュレーションを行っているため、実際の材料状態を正確に捉えることができ、粘弾性材料の非線形性を示す物性値が、材料の変形の速度、大きさによって非線形的に変化する現象を精度良く表現することができる。
【0101】
さらに、シミュレーションに用いられるひずみ、ひずみ速度、応力の値は、粘弾性材料からなる製品の実使用状態を想定した条件下で測定されているため、実際の粘弾性材料の種々の変形状態に対応したシミュレーションを行うことができる。
【0102】
従って、材料の変形状態により、ひずみ、ひずみ速度の関係が変化し、損失係数等の材料物性が非線形性を示すような粘弾性材料からなる製品においても、シミュレーションにより性能や製品の動的挙動を精度良く解析することができる。
【0103】
実際にゴルフクラブヘッドで、ゴルフボールが打撃されたときと同等のひずみ、ひずみ速度の条件下で、実際の打撃試験を正確にシミュレートすることができ、実際の打撃に近い状態でのゴルフボールの性能や変形挙動を精度良く予測することができる。これにより、ゴルフボールの性能の優劣を左右する材料物性の把握が容易となり、製品の性能向上に役立つだけでなく、ゴルフボールの設計段階において、実際のボール試作回数を減らし、試作に要する費用と時間を削減することができる。
【図面の簡単な説明】
【図1】 ひずみεの時刻歴データを示す図である。
【図2】 ひずみ速度ε’の時刻歴データを示す図である。
【図3】 応力σの時刻歴データを示す図である。
【図4】 ひずみ−応力線図と縦弾性係数の算出方法を示す図である。
【図5】 本実施形態の粘弾性モデルとして用いた2要素のVoigtモデルを示す図である。
【図6】 粘性抵抗の時刻歴データを示す図である。
【図7】 ひずみと粘性抵抗の関係を示す図である。
【図8】 ひずみとひずみ速度の関係を示す図である。
【図9】 ゴルフボールモデルのメッシュによる分割状況を示す図である。
【図10】 アルミニウム製中空棒モデルのゴルフボールモデルへの衝突状況を示し、(A)は衝突前、(B)は衝突中、(C)は衝突後の図である。
【図11】 損失係数の算出方法を示す図である。
【図12】 スプリットホプキンソン棒試験機が示された模式的正面図である。
【図13】 試験片のひずみ時刻歴の状態が示された図である。
【符号の説明】
10 ゴルフクラブモデル
11 要素
12 節点
20 アルミニウム製中空棒モデル[0001]
BACKGROUND OF THE INVENTION
The present invention relates to a simulation method for predicting the performance of a product made of a viscoelastic material, and more particularly to a simulation method for accurately predicting the performance of a product made of a viscoelastic material by simulation.
[0002]
[Prior art]
Viscoelastic materials represented by polymer materials such as rubber and elastomer are widely applied to various products such as tires, various balls used in sports competitions, and rolls used in printing presses.
[0003]
For various industrial products such as viscoelastic materials and metal materials as described above, product development using simulation is being conducted in various fields in order to save costs and time for trial production. For example, in order to predict the resilience performance of a golf ball, an actual impact test simulation method has been proposed.
[0004]
Conventionally, when performing the simulation as described above, the physical property values of the constituent materials of the ball measured by a viscoelastic spectrum meter that generally measures the rigidity and viscosity of the material, a tensile tester that measures the longitudinal elastic modulus (Young's modulus), etc. Is used as an input value to the simulation. In particular, the viscoelastic spectrum meter is useful for simulating a product made of a viscoelastic material because a physical property value of a test piece to which dynamic strain is applied is measured.
[0005]
[Problems to be solved by the invention]
However, the measurement with the above-mentioned viscoelastic spectrum meter, tensile tester for measuring the longitudinal elastic modulus, etc. cannot give a large amount of deformation to the test piece, and the maximum strain rate applied to the test piece made of viscoelastic material at the time of measurement. Is as small as about 0.001 / s to about 1.0 / s, and the maximum compressive strain is as small as about 0.0001 to 0.02.
[0006]
On the other hand, a product made of a viscoelastic material may cause a large deformation at high speed due to the influence of external force during actual use. For example, in the material constituting the golf ball, the maximum strain rate is about 500 / s to 5000 / s and the maximum compressive strain is mainly about 0.05 to 0.50 at the time of actual hitting. Large value.
[0007]
As described above, in the measurement using a viscoelastic spectrum meter, a tensile tester for measuring the longitudinal elastic modulus, the physical property value of the viscoelastic material is measured under the same high-speed and large deformation condition as the actual use condition as described above. The maximum strain rate and the maximum compressive strain are greatly different from those in actual use. For this reason, in the conventional simulation which inputs the material physical property value obtained by the viscoelasticity spectrum meter, the tensile testing machine, etc., there exists a problem that the accurate simulation which considered the physical property of the viscoelastic material cannot be performed.
[0008]
That is, it is known that the deformation behavior of a viscoelastic material subjected to an impact load is significantly influenced by the deformation amount or the deformation speed, unlike the case of receiving a static load. In particular, a polymer material typified by rubber or elastomer causes deformation at a very high speed of 10,000 minutes or several thousandths of a second due to an impact load, and the deformation amount is as large as several tens of percent. There are many viscoelastic materials with such high-speed and large-deformation behavior, and high-precision simulation is required for efficient product development. Specifically, in products with impact during use, such as golf balls, the dynamic behavior and material characteristics under high-speed and large-deformation conditions affect the performance of the product, so accuracy under actual use conditions is good. Simulation is essential for product development.
[0009]
In addition, viscoelastic materials include materials whose physical properties such as loss factors change depending on the magnitude of strain and strain rate when an external force such as impact load is applied. That is, since the viscoelastic material has a wide range of deformation speed and size, the physical property value is not linear but has a strong non-linear change depending on the deformation speed and size. Specifically, as the material is deformed by an external force and the strain increases, the loop area of the SS (strain-stress) curve increases, and the physical property value such as a loss coefficient is determined by the deformation state (deformation of the material). It changes with speed, size), etc. and shows non-linearity. Depending on the viscoelastic material, there are many materials having strong non-linearity, and a highly accurate simulation is required for a product using such a viscoelastic material.
[0010]
However, there has been no method capable of accurately expressing a physical property indicating nonlinearity of a viscoelastic material, for example, a phenomenon in which a loss factor strongly and nonlinearly changes depending on the deformation speed and size of the material. Conventionally, for products such as golf balls made of viscoelastic materials, simulation has been performed on the assumption that viscoelastic property values hardly change. As a result, there is a problem that the performance of a product made of a viscoelastic material at the time of actual use cannot be accurately predicted by simulation. In fact, the performance of the product must be evaluated by producing a prototype. There is no problem.
[0011]
The present invention has been made in view of the above-described problems, and it is an object of the present invention to accurately predict the performance under actual use conditions of a product made of a viscoelastic material whose material properties exhibit nonlinearity.
[0012]
[Means for Solving the Problems]
In order to solve the above-mentioned problems, the present invention measures the strain, strain rate, and momentary values of stress generated in the viscoelastic material under measurement conditions assuming the actual use state of a product made of a viscoelastic material,
Based on the time history data of the strain, strain rate, and stress, and the viscoelastic model taking into account the viscosity of the viscoelastic material, the viscoelasticity model is divided into cases when the strain is loaded and when the strain is unloaded or restored. Deriving time history data of viscous resistance of elastic materials,
The product made of the viscoelastic material is set as a product model to be analyzed, the product model is divided into a number of elements, the relationship between the strain, strain rate, and viscous resistance is input, and when the strain is loaded When analyzing the finite element method in consideration of the difference in viscous resistance between unloading or restoring the strain and predicting the performance of the product model made of the above viscoelastic material,
The strain and strain rate of the global coordinate system generated for each element are decomposed into a deviation component and a volume component, and the strain and strain rate of the deviation component are converted from the global coordinate system to the main strain coordinate system and the main strain rate coordinate system. ,
For each coordinate axis in the main strain coordinate system and the main strain rate coordinate system using the converted deviation main strain and deviation main strain rate, even when the strain is the same, the strain is loaded and the strain is unloaded or restored. A simulation method for predicting the performance of a product made of a viscoelastic material is provided.
[0013]
In this way, the strain and strain rate measured under actual use conditions are decomposed into volume components and deviation components, and converted to the main strain coordinate system and the main strain rate coordinate system. For each axis of the main strain coordinate system and the main strain rate coordinate system, even when the strain is the same value, a different viscous resistance is determined depending on whether the strain is loaded or unloaded or restored. The resistance is variable for each element in three axial directions for each element. For this reason, whether the deformation state of the material is a load state or an unloading state, it is possible to perform an accurate analysis corresponding to the actual situation, and to perform an accurate analysis in the three-dimensional direction, The performance prediction accuracy of simulation can be improved.
[0014]
That is, the deformation state occurring in the viscoelastic material can be divided into a “load state” in which the applied strain increases in the compression direction and a “restored state” in which the amount of compression gradually decreases. The strain rate relationship is different. Therefore, as described above, the strain state generated in the viscoelastic material takes into consideration two states, the load state and the unloading (or restoration) state. It is possible to accurately determine the actual deformation state of the material by determining the viscous resistance using the relationship and deriving different viscous resistance values in the strain loading state and unloading state even when the strain is the same value. it can. Therefore, it is possible to accurately represent the phenomenon that the material changes depending on the deformation speed, the size, and the like, and the simulation accuracy can be further improved.
[0015]
Also, by calculating the magnitude (norm) of the principal strain for each element and comparing it with the norm of the previous step at the time of simulation, if the norm increases, the load, if the norm decreases It is preferable to determine unloading and determine the strain load and the unloading state.
For example, with only one direction of main strain, the load and unloading show various changes with a small change. Thus, by making the size in multiple directions, the load and unload must be changed to some extent. Since unloading does not change, calculation stability can be improved.
[0016]
In the present invention, ε in the global coordinate systemx, Εy, Εz, Εxy, Εyz, ΕzxIn order to consider viscoelastic characteristics in the deviation component of strain and strain rate generated for each element, the above strain and strain rate are decomposed into volume component and deviation component, and converted to main strain coordinate system and main strain rate coordinate system Using the principal strain and principal strain rate of the measured deviation component, the viscous resistance that is different between the strain loading state and the strain unloading or restoring state is determined for each axis of the principal strain coordinate system and the principal strain rate coordinate system. The viscous resistance is variable for each element in each of the three axial components. For this reason, the analysis which considered the viscoelasticity of the viscoelastic material correctly can be performed about a three-dimensional direction, and the precision of a simulation can be improved.
Here, the volume component is εv= (Εx+ Εy+ Εz), And the deviation component is εxy′ = Εxy/ 2, εyz′ = Εyz/ 2, εzx′ = Εzx/ 2, εx′ = Εx−εv/ 3, εy′ = Εy−εv/ 3, εz′ = Εz−εv/ 3 refers to the component represented. The reason why the shear component is halved is that the physical strain is converted.
The reason why only the deviation component is converted is that the viscoelastic component of the volume component is very small and cannot be obtained by actual measurement at present. That is, since there is almost no viscoelastic component of the volume component, judgment is made only by the deviation component, and since the transformation has only a result in one axis direction from the actual measurement result, εxy, Εyz, ΕzxConvert so that component = 0, ε1, Ε2, Ε3The actual measurement results are used for the three components (main strain three components) only.
The global coordinate system refers to the initial shape, that is, the coordinate system that determines the shape of the model, and the main strain coordinate system and the main strain rate coordinate system have strain and strain rate with a shear component of zero. When converted in this way, it indicates a coordinate system obtained by converting the entire coordinate system by the same angle.
[0017]
Although the strain (strain rate) generated for each element is 6 components, the material physical property measured under the condition of high speed and large deformation is the material physical property in the uniaxial direction. Is insufficient. Accordingly, the six strain components generated for each element are converted into the main strain coordinate system, and the three strain components in the main axis direction are obtained.
That is, by reducing the strain in the shear direction to 0, the six components (εx, Εy, Εz, Εxy, Εyz, Εzx) Is the strain in the principal axis direction (principal strain coordinate system, ε1, Ε2, Ε3) 3 components. Thereby, the material characteristic of a viscoelastic material can be evaluated from an actual measurement result, without considering a shear direction.
[0018]
Specifically, the strain and strain rate of the entire coordinate system generated for each element are decomposed into a deviation component and a volume component, and viscoelasticity is considered in the deviation component. Dissolved strain and strain rate deviation 6 components are converted from the global coordinate system to the main strain coordinate system and the main strain rate coordinate system, and the deviation main strain and the main strain rate deviation in the main strain coordinate system and the main strain rate coordinate system are converted. These three components are used. Thereby, viscous resistance can be determined for each axis in the main strain coordinate system and the main strain rate coordinate system.
[0019]
In addition, considering the viscoelasticity from the above viscous resistance, the stress and strain of the deviation component are calculated, the stress and strain of the deviation component are reconverted to the original global coordinate system, and the stress and strain of the volume component and the deviation component are calculated. The stress and strain for each element are calculated from the stress and strain. Thus, by calculating the stress and strain for each element, it is possible to calculate the stress and strain every moment.
[0020]
Specifically, using the deviation principal strain and the deviation principal strain rate, the different viscous resistances are determined for each coordinate axis when the strain is loaded and when the strain is unloaded or restored. Then, the stress and strain of the deviation component considering viscoelasticity are calculated from the viscous resistance.
(Formula 1)
Figure 0003910825
Where Ei: Deviation principal strain of i-th element, n: analysis time step,
βi: Β of the i-th element (where β = Et/ Η, Et= Eshort-Elong, Η: viscous resistance), γi: Indicates the deviation principal strain rate of the i-th element.
(Formula 2)
Figure 0003910825
Where σi: I-th element stress, Fi: Rigidity coefficient, Gshort: Short-term shear modulus, Glong: Indicates long-term shear coefficient.
After calculating β and stiffness coefficient F in Formulas 1 and 2 from the experimental data, the strain deviation 6 component is updated by Formula 1, and the stress is calculated by Formula 2. Thereafter, the principal strain coordinate system and the main strain rate coordinate system are converted to the entire coordinate system, and the stress is updated by the following Equation 3.
(Formula 3)
Figure 0003910825
Where K: bulk modulus, εv: Volumetric strain, σij: Refers to the stress in the global coordinate system.
[0021]
In addition, using a viscoelastic model that takes into account the viscosity of the viscoelastic material, a different viscous resistance is derived between the strain loading state and the strain unloading or restoring state, and the strain, strain rate, and strain loading state are derived. Since the simulation is performed by inputting into the product model the relationship between the viscous resistance that differs depending on whether the strain is unloaded or restored, the physical property value indicating the nonlinearity of the viscoelastic material is nonlinear depending on the deformation speed and magnitude of the material. Phenomenon that changes automatically can be expressed with high accuracy. In addition, the values of strain, strain rate, and stress used in the simulation are measured under conditions that assume the actual usage state of products made of viscoelastic materials, so they correspond to various deformation states of actual viscoelastic materials. Simulation can be performed. Therefore, even for products made of viscoelastic materials in which the relationship between strain and strain rate changes depending on the deformation state of the material, and the material properties such as loss coefficient show nonlinearity, it is possible to accurately predict the performance by simulation. ing.
[0022]
In this simulation method, the values of strain, strain rate, and stress that occur in the viscoelastic material are measured every moment under measurement conditions that assume the actual usage state of the product made of the viscoelastic material. Specifically, the measurement conditions are determined assuming that an external force is applied to the product during actual product use and the viscoelastic material is deformed. Under the above measurement conditions, strain, strain rate, and stress values generated in the viscoelastic material are measured every moment to obtain each time history data. Therefore, information on the deformation state of the viscoelastic material when an external force assuming an actual use state is applied to the viscoelastic material can be obtained from the time history data of the strain, strain rate, and stress. Thereby, the material physical property can be accurately predicted even for a viscoelastic material having a high-speed large deformation behavior due to an impact load.
[0023]
It is preferable to measure the values of the strain, strain rate, and stress every moment under a plurality of different measurement conditions. By changing the magnitude of the external force applied to the product and setting the measurement conditions for multiple patterns, it is possible to obtain various pattern data for strain, strain rate, and stress, and improve the accuracy of simulation input values. be able to. In addition, in order to obtain as much pattern data as possible, the strain, strain rate, and stress values are measured until external strain is applied to the viscoelastic material and the strain is almost zero. It is preferable that the measurement time interval is shorter than the simulation accuracy.
[0024]
As the viscoelastic model considering the viscosity of the viscoelastic material, a viscoelastic model including a spring and a dashpot is preferable. Such a viscoelastic model can simplify the viscosity of the viscoelastic material, so that the influence of the viscosity on the deformation state of the material can be easily taken into consideration. Specifically, a Maxwell model, a Vogt model, or a combination of a plurality of springs and dashpots can be used. A two-element model is preferable from the viewpoint of simplifying the model. In such a viscoelastic material model, the viscosity resistance η of the dashpot and the stiffness of the spring (longitudinal elastic coefficient E or shear coefficient (lateral elastic coefficient)) are made variable. The shear coefficient is a physical property value determined by a longitudinal elastic modulus (Young's modulus) E and Poisson's ratio.
[0025]
Based on the above time history data of strain, strain rate, and stress and viscoelastic model considering the viscosity of the viscoelastic material, the viscoelastic material has different viscosities when the strain is loaded and when the strain is unloaded or restored. Resistance time history data is derived. Specifically, according to the above viscoelastic model, the strain, strain rate, and longitudinal elastic modulus of the viscoelastic material, the viscous resistance that is different between the strain loading state and the strain unloading or restoring state. Each of these relationships is formulated. By formulating in this way, viscous resistance is expressed as a function of strain and strain rate. Then, the longitudinal elastic modulus of the viscoelastic material is derived from the strain and stress generated in the viscoelastic material obtained by the above measurement, and substituted into the above function together with the strain and strain rate, so that the strain load condition and the strain can be removed. Different values of viscous resistance are derived for the loaded or restored state. Since time history data on strain, strain rate, and stress is obtained by the above measurement, corresponding time history data can be obtained similarly for viscous resistance. Thus, the viscous resistance is determined by the values of strain and strain rate, and the viscous resistance also changes corresponding to the time change.
[0026]
A product made of the viscoelastic material is set as a product model to be analyzed, and the strain, strain rate, and strain load are included in the input data for calculation (or input data) including the product model and including speed, constraint conditions, etc. A simulation is performed in consideration of the change in viscous resistance due to the difference in strain and strain rate by inputting the relationship between the viscous resistance that is different between the state and the state of unloading or restoring the strain.
In the simulation calculation, the relationship between the above-mentioned strain, strain rate, and the viscous resistance that is different between the strain loading state and the strain unloading or restoring state is entered in the product model. And the two-dimensional data of the relationship between strain and viscous resistance can be input and calculated. Also, the relationship between strain, strain rate, and viscous resistance can be input as three-dimensional data, and the viscous resistance can be input as a function and calculated.
[0027]
The two-dimensional data of the relationship between strain and strain rate, and the relationship between strain and viscous resistance are obtained by using the above-mentioned relationships, when strain, strain rate and corresponding strain are loaded, and strain is unloaded or restored. Different viscous resistances are written as input data. Specifically, strain, strain rate, etc. are measured under multiple measurement conditions, and for each pattern with different measurement conditions, the corresponding relationship between strain and strain rate is recorded from the time series data of strain and strain rate. To do. Also record the value of viscous resistance corresponding to each curve. By arranging these three correspondences of strain, strain rate, and viscous resistance, the viscous resistance at a certain strain and strain rate is accurately derived and calculated.
[0028]
As the number of measurements under different strain and strain rate conditions increases, the physical properties of the viscoelastic material can be realized with higher accuracy, so it is possible to measure strain, strain rate, etc. under multiple different measurement conditions as described above. preferable. However, the more measurement data, the more time required for simulation, so if the strain and strain rate are not the same as the measurement data measured under the specified measurement conditions, the viscous resistance is calculated using interpolation. Is preferred. The above-mentioned interpolation can be performed by various methods such as linear interpolation from two values of viscous resistance under conditions of close strain and strain rate, or interpolation using individual measured values under all measured conditions. Can be performed. By performing such an interpolation operation, it is possible to calculate the viscous resistance corresponding to changes in strain and strain rate generated in the material due to differences in measurement conditions.
[0029]
Even if the viscoelastic material is a material whose material properties exhibit nonlinearity, according to the present simulation method, the material properties and deformation behavior assuming an actual use state can be accurately simulated. Even in a material in which the material property of the viscoelastic material is not linear, and exhibits strong nonlinearity that varies depending on the deformation speed and the magnitude of the deformation, the viscous resistance is distorted as described above. Nonlinearity can be taken into account by calculating the viscous resistance using a viscoelastic model determined by two values. In this simulation method, it is possible to predict the performance of a product assuming an actual use state by performing a simulation with high accuracy, particularly in a material having nonlinearity with a strong loss coefficient.
[0030]
The simulation method of the present invention is performed by the finite element method as described below.
When the simulation is performed by finite element analysis, a large number of nodes and elements are set in the product model. That is, when predicting material properties by simulating a product composed of viscoelastic materials using the finite element method, the stiffness (longitudinal elastic modulus or shear coefficient) of the viscoelastic model spring is determined, and viscous resistance indicating viscosity As described above, by determining the strain and strain rate generated in each element for each element, the material physical properties under conditions of proper strain and strain rate can be expressed for each element. Instead of inputting the longitudinal elastic modulus, a shear coefficient may be input from the relationship with the Poisson's ratio. It can be changed according to the specification of the finite element method program whether to input the longitudinal elastic modulus or the shear modulus.
[0031]
The strain, strain rate, and stress are measured by a split Hopkinson bar tester. Using a split Hopkinson bar tester, the specimen can be strained at a high speed and with a large deformation, at a high speed of 10,000 minutes or a few thousandths of a second, and a deformation amount of several tens of percent. Time series data of strain, strain rate, and stress of viscoelastic material can be obtained. If the measurement conditions of the split Hopkinson bar tester are the same conditions as the strain and strain rate generated in the viscoelastic material when an impact load is applied to the product, the high-speed large deformation state of the viscoelastic material caused by external force, etc. The material physical properties of the viscoelastic material corresponding to various states can be obtained. Therefore, the accuracy of the simulation can be improved by using the material properties measured by the split Hopkinson bar tester.
[0032]
In addition, the split Hopkinson bar testing machine can easily measure the material physical properties of the specimen in various strain and strain rate areas simply by changing the impact speed of the striking bar that applies impact to the specimen. In addition, material properties can be obtained in various strain and strain rate patterns.
[0033]
The split Hopkinson bar tester was originally used for the evaluation of impact behavior of metallic materials. In the present invention, the split Hopkinson bar tester is used to improve the evaluation of viscous viscoelastic materials. ing. The measuring method of the split Hopkinson bar tester will be described later.
[0034]
If the specimen can be strained at a high speed and with large deformation, and the physical properties of the material can be measured under the conditions of strain and strain rate under the conditions in which the product is actually used, the above split Hopkinson can be used. It goes without saying that the material properties may be obtained by a measuring method other than a bar.
[0035]
The maximum value of the strain generated in the viscoelastic material when measuring strain, strain rate, and stress under the measurement conditions assuming the actual use state of the product made of the viscoelastic material is in the range of 0.05 to 0.50. Or / and the maximum value of the strain rate is in the range of 500 / s to 10000 / s, preferably 500 / s to 5000 / s. The range between the maximum value of strain and the maximum value of strain rate is the condition of strain and strain rate generated during high-speed large deformation of viscoelastic material. It is preferable to use three time series data of strain, strain rate and stress under the conditions.
[0036]
The viscoelastic material is a golf ball material, and the product model is a golf ball. A golf ball is made of a viscoelastic material and is subjected to external force such as impact load during actual use and is subject to high-speed deformation or large deformation. The high-speed and large deformation state greatly affects the performance of the golf ball itself. Effect. For this reason, the analysis by the simulation method of the present invention is very useful for predicting the performance of the golf ball, and can accurately predict the performance of the golf ball without making a prototype. Can be achieved.
[0037]
A collision phenomenon between the golf ball and a hitting object assuming a golf club head is simulated, and the behavior of the golf ball at the time of the collision is predicted. Since the golf club head can predict the material properties of the viscoelastic material constituting the golf ball under the same strain, strain rate, and stress conditions as when the golf ball was hit, It is possible to predict the behavior of the golf ball, such as the coefficient of restitution and the deformation state of the golf ball when the golf ball is hit.
[0038]
The golf ball may be a so-called one-piece golf ball made of a single substance such as a crosslinked rubber layer, or may be a so-called two-piece golf ball having a core covered with a cover, such as a crosslinked rubber layer, and more than three layers. A so-called multi-piece golf ball may be used, and can be applied to golf balls having any structure made of a viscoelastic material.
[0039]
The viscoelastic material includes all materials having viscoelasticity. For example, a thermoplastic resin, a thermosetting resin, various elastomers, rubber, etc. are mentioned, These single bodies or a mixture can be used. Moreover, various additives, such as a coloring agent, a deterioration inhibiting agent, a crosslinking agent, may be mix | blended with these simple substance and mixture as needed. Any material that can measure strain, strain rate, and stress under the actual usage conditions of the product can be applied.
[0040]
Specifically, synthetic resins such as ionomer resins used as a constituent material of golf balls, polybutadiene (butadiene rubber), natural rubber, polyisoprene, styrene-butadiene copolymer, ethylene-propylene-diene copolymer (EPDM) And urethane rubber.
[0041]
Examples of the product made of the viscoelastic material include a rubber roller, a rubber belt, or a tire used in a printing machine such as a printer other than a golf ball, and sports equipment for tennis and golf made of other viscoelastic materials. As a product made of viscoelastic material, it is sufficient that the viscoelastic material is used for at least a part of the product, and it may be a composite molded product with other materials such as metal materials. It is also possible to predict the performance of the part corresponding to. In particular, it is suitable for a product that receives an impact load and exhibits a large amount of deformation at high speed during use, and the performance and dynamic behavior of the product can be accurately predicted by this simulation method.
[0042]
In this simulation method, the strain, strain rate, longitudinal elastic modulus or shear coefficient is input to the calculation input data including the product model, and the longitudinal elastic modulus or shear coefficient changes due to the difference in strain and strain rate. By taking into account, viscoelastic materials whose longitudinal elastic modulus varies depending on the deformation state of the material (strain, strain rate, etc.) can be analyzed accurately, and simulation accuracy can be further improved. it can.
Moreover, depending on the case, it can also be classified according to the stress loading state and the unloading state.
[0043]
DETAILED DESCRIPTION OF THE INVENTION
Hereinafter, embodiments of the present invention will be described with reference to the drawings.
In the simulation method of the first embodiment of the present invention, a material mainly composed of butadiene rubber, which is a constituent material of a golf ball, is used as a viscoelastic material exhibiting nonlinearity. By using a split hopkinson bar testing machine improved for measuring viscoelastic materials, using the butadiene rubber as a main component as a test piece, the measurement conditions at the time of large deformation at high speed assuming the actual use state of the golf ball, The strain, strain rate, and stress values that occur in materials based on butadiene rubber are measured. A measuring method using a split Hopkinson bar tester will be described later.
[0044]
The split Hopkinson bar tester can measure material properties in various strain and strain rate regions by changing the impact speed of the striking bar. In this embodiment, four patterns of collision speeds (7 m / s, 14 m / s, 20 m / s, 25 m / s) are employed, and measurement is performed under four measurement conditions. Time history data of speed and stress is obtained. The time history data of the strain ε, strain rate ε ′, and stress σ are shown in FIG. 1, FIG. 2, and FIG. 3, respectively.
[0045]
After striking the striking rod with the split Hopkinson bar testing machine, as shown in FIG. 1, the material mainly composed of butadiene rubber is strained, and the strain value increases until a certain time, It is gradually decreasing. As shown in FIG. 2, the strain rate shows a positive value until the strain reaches the maximum value, and shows a negative value after the maximum strain value. As shown in FIG. 3, the value of stress similarly changes with time.
[0046]
A strain-stress diagram is drawn from the strain time history data and stress time history data as shown in FIG. In the strain-stress diagram of FIG. 4, the strain ε when the strain value is the maximummaxAnd the value of the stress σa at that time, the longitudinal elastic modulus E of the test piece is calculated for each collision pattern by the following mathematical formula 4.
(Formula 4)
E = σa/ Εmax      ---- (4)
In addition, in the material mainly composed of butadiene rubber used in this embodiment, the longitudinal elastic modulus in each collision pattern shows almost the same value even if the impact speed of the hitting bar is changed in the split Hopkinson bar tester. Yes.
[0047]
For products made of viscoelastic materials, a viscoelastic model that takes into account the viscosity of the viscoelastic material is set in order to perform a simulation in consideration of the viscosity. Specifically, in this embodiment, a basic two-element Voigt model is used as a viscoelastic model including a spring and a dashpot as shown in FIG. That is, the viscoelastic material model uses the dashpot viscous resistance η and the stiffness of the spring (longitudinal elastic modulus E or shear coefficient) as variable.
[0048]
As shown in FIG. 5, in the viscoelastic model that is a two-element Voigt model,1, Σ is the stress generated in the dashpot2Then, the stress σ generated as a whole can be expressed by the following Equation 5.
(Formula 5)
σ = σ1+ Σ2
= Eε + ηε '--- (5)
[0049]
Therefore, from the above formula 5, the viscous resistance η can be expressed by the following formula 6 in such a viscoelastic model.
(Formula 6)
η = (σ−Eε) / ε ′ −−− (6)
[0050]
Accordingly, the viscous resistance η corresponding to the strain ε and the strain rate ε ′ is calculated from time to time using the time history data of the strain ε, strain rate ε ′, and stress σ shown in FIGS. it can. Here, as described above, the values of the strain, strain rate, and stress change with the passage of time, and the strain state generated in the viscoelastic material is the “load state in which the applied strain increases in the compression direction. "And a" restoration state "in which the amount of compression gradually decreases. Therefore, when calculating the viscous resistance, the strain state generated in the viscoelastic material is calculated separately in two states: a load state and an unloading (or restoration) state, as shown in FIG. The time history data of viscous resistance is obtained at each collision speed.
[0051]
Thus, after determining the longitudinal elastic modulus E of the viscoelastic material from the strain and stress obtained by the measurement under the high speed and large deformation condition, the above formula 6 derived from the viscoelastic model taking viscosity into consideration, Viscosity resistance η is calculated separately during loading and unloading (or restoring) states, and different viscous resistances η are applied to the simulation when loading and unloading or restoring strain. To do.
[0052]
FIG. 7 shows the relationship between strain and viscous resistance based on the strain time history data and viscous resistance time history data. As shown in FIG. 7, the value of the viscous resistance shows a different value between the strain loading state and the unloading (or restoration) state even when the strain is the same value. FIG. 8 shows the relationship between strain and strain rate based on strain time history data and strain rate time history data. Thus, the viscous resistance changes depending on the strain and strain rate, and shows different values in the strain loading state and unloading (or restoration) state. Therefore, the viscous resistance is determined by the strain and strain rate. The viscous resistance is a function of strain and strain rate, which is divided into cases of strain loading and unloading. In FIG. 8, the data of the four collision velocities appear to overlap in a straight line on the display in the negative strain velocity area, but the actual values are different for each collision velocity. It is the speed value.
[0053]
The relationship between the strain and viscous resistance shown in FIG. 7 and the relationship between the strain and strain rate shown in FIG. 8 are written as input data in each collision pattern as will be described later, and the strain and strain are analyzed by finite element analysis. A simulation is performed in consideration of a change in viscous resistance due to a difference in speed, that is, a change in viscous resistance due to a difference between a strain loading state and an unloading (or restoration) state.
[0054]
In the present embodiment, a golf ball model is set as a product model made of a viscoelastic material, and simulation is performed assuming a collision of a golf club head (striking object) with a golf ball. As shown in FIG. 9, the golf ball model 10 assumes a one-piece ball mainly composed of butadiene rubber and has a spherical shape with a diameter of 42.8 mm.
[0055]
Initial conditions are set in the product model for finite element analysis. That is, initial conditions such as the size, shape, and structure of the golf ball model 10 are set, and the golf ball model 10 is divided into a large number of elements 11 to obtain a large number of nodes 12. The total number of elements in the entire range of the golf ball model is preferably 1000 to 100,000, and more preferably 2500 to 20000. Note that the upper limit value is based on the ability of the computer at the present stage, and as the ability of the computer improves in the future, the analysis time will be shortened, so it is easy to imagine that the upper limit range will change. it can.
[0056]
Based on the set conditions, the relationship between the strain, strain rate, and viscous resistance and the data of the golf ball model 10 are written in the input data during the simulation calculation. When the calculation is executed, an appropriate viscous resistance is calculated for each element from time to time, and the formula is calculated using the viscous resistance. In this simulation, two-dimensional data of the relationship between strain and strain rate and the relationship between strain and viscous resistance divided according to the strain loading state and unloading (or restoration) state are input and calculated. Yes. From the time series data of strain and strain rate for each of the four patterns with different measurement conditions, record the corresponding relationship between strain and strain rate, and write the strain, strain rate and corresponding viscous resistance as input data. Become. In the case of strains and strain rates that are not the same as the strain and strain rate measurement data measured under the four measurement conditions, linear interpolation is performed using two values of the viscous resistance under the conditions of adjacent strains and strain rates. ing.
[0057]
More specifically, when one element is focused on, the norm is increased by calculating the magnitude (norm) of the strain generated in that element and comparing it with the norm of the previous step during simulation. If the load and norm are reduced, it is determined that the load is unloaded, the strain load and the unload state are determined, and information on the strain and strain rate in consideration of the material state is obtained. Next, from the relationship between the strain and strain rate of each collision case (measurement case) obtained from the measurement, refer to the value of the strain rate when the same strain in each case occurs, and Search for two cases with strain rate values close to the strain rate value. By using the strain rate and viscous resistance values when the strains corresponding to these two cases are the same value, it is possible to calculate the appropriate viscous resistance value corresponding to the strain and strain rate generated by the element by interpolation. it can. In this embodiment, interpolation by simple primary interpolation is performed in this way. In this way, by determining the viscous resistance for loading and unloading depending on the strain and strain rate generated for each element, the material physical properties under the conditions of appropriate strain and strain rate for each element. Represents.
[0058]
Also, the strain and strain rate of the global coordinate system generated for each element are decomposed into a deviation component and a volume component, and the strain and strain rate of the deviation component are changed from the global coordinate system to the main strain coordinate system and the main strain rate coordinate system. Using the converted deviation main strain and deviation main strain rate, as described above, for each coordinate axis, as shown above, even when the strain is the same, the strain is loaded and the strain is unloaded or restored. Viscosity resistances different from each other are determined, and stresses and strains of deviation components taking viscoelasticity into account are calculated from the viscous resistances by the formulas 1, 2, and 3. Thereafter, the stress and strain of the deviation component are reconverted into the entire coordinate system, and the stress and strain for each element are calculated from the stress and strain of the volume component and the stress and strain of the deviation component.
[0059]
That is, the six strain components generated for each element are converted into the main strain coordinate system, and the three strain components in the principal axis direction are obtained. The distortion and strain rate of the entire coordinate system generated for each element are decomposed into a deviation component and a volume component, and viscoelasticity is taken into consideration for the deviation component. Dissolved strain and strain rate deviation 6 components are converted from the main coordinate system to the main strain coordinate system and the main strain rate coordinate system, and the deviation main strain and the main strain rate deviation in the main strain coordinate system and the main strain rate coordinate system are converted. Using each of the three components, a different viscous resistance is determined for each axis in the main strain coordinate system and the main strain rate coordinate system, even when the strain is the same, depending on whether the strain is loaded or unloaded or restored. is doing.
[0060]
In the present embodiment, specifically, as shown in FIGS. 10A, 10B, and 10C, the golf ball model 10 is made of a cylindrical aluminum having a hitting object of 200 g (same as the weight of the golf club head). The material properties of the golf ball model 10 when the hollow rod model 20 is collided at a speed of 45 m / s are analyzed by a finite element method and simulation is performed. As a result, the amount of strain generated in each element 11 of the golf ball model 10 during a predetermined time from when the hit object collides to after a predetermined time is calculated.
[0061]
Mass is distributed to each node constituting one element, each node is replaced with a mass point, and the speed of each node is divided as the speed of the mass point to obtain the total speed. That is, the velocity Vbi (i = x, y, z) of the ball after impact is calculated as follows using Equation 7. The total momentum of the ball is considered to be the sum of the momentum of all mass points, and the total momentum divided by the total weight is defined as the ball jumping speed Vbi.
[0062]
(Formula 7)
Figure 0003910825
Here, N and M are the number of nodes of the ball, the total mass, and Vni and Mn are the nth translational speed and the mass of the element including the node divided by the number of nodes included in the element.
[0063]
The circular surface 20a of the aluminum hollow rod model 20 is a collision surface, the collision surface is flat, and the golf ball model 10 is a flat frontal collision. The center point 20 b of the circular surface 20 a of the aluminum hollow rod model 20 is caused to collide with the first golf ball model 10.
[0064]
By the method described above, the speeds of the aluminum hollow rod model 20 and the golf ball 10 before and after the collision are calculated, and the coefficient of restitution of the golf ball 10 is calculated and predicted from the respective speeds and weights.
[0065]
As described above, the strain and strain rate are decomposed into a volume component and a deviation component, and converted into the main strain coordinate system and the main strain rate coordinate system, and the main strain coordinate system and the main strain rate are used. For each axis of the strain rate coordinate system, even when the strain is the same, a different viscous resistance is determined in the strain loading state and in the strain unloading or restoring state, and the viscous resistance is determined in each axial direction 3 for each element. The components are variable. For this reason, the analysis which considered the viscoelasticity of the viscoelastic material correctly can be performed about a three-dimensional direction, and the precision of a simulation can be improved.
[0066]
Also, the viscous resistance divided between the strain loading state and the unloading (or restoration) state is determined as a function of strain and strain rate, and the relationship between the above-mentioned strain, strain rate, and viscous resistance is golfed. It is input to the ball model and simulated by finite element analysis. For this reason, the corresponding viscous resistance can be calculated from the momentary strain and strain rate of each element, and the viscous resistance value takes into account the difference between the strain loaded state and unloaded (or restored) state. The properties of viscoelastic materials having strong nonlinearity can be predicted with high accuracy even under high speed and large deformation conditions.
[0067]
Therefore, it is possible to accurately grasp the performance, dynamic behavior, etc. of a golf ball model made of a viscoelastic material under the same conditions as when hit with an actual golf club head, without making a prototype. .
[0068]
Further, as shown in FIG. 11, the phase angle δ is calculated from the stress-strain curve using the following formula 8, and the loss factor (tan δ) can also be calculated from the phase angle δ.
(Formula 8)
δ = sin-1((Σa−σb) / Σmax---- (8)
[0069]
(Measurement of material properties by split Hopkinson bar tester)
FIG. 12 is a schematic front view showing a split Hopkinson bar testing machine improved for measuring viscoelastic materials.
[0070]
The split Hopkinson bar testing machine shown in FIG. 12 includes a striking bar 51, an input bar 53, and an output bar 55, which are arranged on a straight line. A first strain gauge 57 and a second strain gauge 59 are attached to the input bar 53. A third strain gauge 61 and a fourth strain gauge 63 are attached to the output bar 55. A cylindrical test piece 70 is sandwiched between the rear end 53a of the input bar 53 and the front end 55a of the output bar 55.
[0071]
The test piece 70 may be formed by molding the viscoelastic material to be measured into the shape of the test piece, or may be cut out from a product state in which the viscoelastic material is molded. In the present embodiment, the length of the test piece 70 (that is, the distance between the rear end 53a of the input bar 53 and the front end 55a of the output bar 55) is 4 mm, and the cross-sectional diameter of the test piece 70 is 18 mm.
[0072]
The striking bar 51, the input bar 53, and the output bar 55 are polymethylmethacrylate cylinders having a cross-sectional diameter of 20 mm, a longitudinal elastic modulus of 5300 MPa, and a specific gravity of 1.19. The length of the hitting stick 51 is 100 mm. The lengths of the input bar 53 and the output bar 55 (hereinafter, the input bar 53 and the output bar 55 are also referred to as “stress bars”) are 2000 mm.
The first strain gauge 57 is attached at a position 900 mm from the rear end 53 a of the input bar 53, and the second strain gauge 59 is attached at a position 600 mm from the rear end 53 a of the input bar 53. The third strain gauge 61 is attached at a position 300 mm from the front end 55 a of the output bar 55, and the fourth strain gauge 63 is attached at a position 600 mm from the front end 55 a of the output bar 55.
[0073]
Thus, the split Hopkinson bar testing machine of FIG. 12 is made of a synthetic resin in which the hitting bar 51, the input bar 53 and the output bar 55 are made of polymethyl methacrylate, and the hitting bar 51 and the input bar 53 are as large as 2000 mm, In addition, since the distance between the first strain gauge 57 and the rear end 53a of the input bar 53 and the distance between the second strain gauge 59 and the rear end 53a of the input bar 53 are large, the viscosity of the crosslinked rubber or the like used for golf balls is large. Suitable for measuring strain, strain rate, and stress of elastic materials.
[0074]
A uniaxial plastic strain gauge is used as the first strain gauge 57, the second strain gauge 59, the third strain gauge 61, and the fourth strain gauge 63. In this embodiment, KFP-5 manufactured by Kyowa Denki Co., Ltd. is used. 350-C1-65 is used and attached to the positions of the input bar 53 and the output bar 55 described above. The attachment positions of the first strain gauge 57 to the fourth strain gauge 63 to the input bar 53 and the output bar 57 are on the same line in the length direction.
[0075]
In the measurement of strain, strain rate, and stress by the split Hopkinson bar testing machine, first, the striking bar 51 collides with the front end 53a of the input bar 53 at a predetermined speed. As a result, an incident distorted wave is generated, and this incident distorted wave travels toward the rear end 53 a of the input bar 53. A part of this incident distorted wave is reflected at the rear end 53 a of the input bar 53, becomes a reflected distorted wave, and travels toward the front end 53 b of the input bar 53. A part of the incident distorted wave passes through the test piece 70 from the rear end 53a of the input bar 53, further propagates to the output bar 55 to become a transmitted distorted wave, and proceeds toward the rear end 55b of the output bar 55.
[0076]
The incident strain wave is actually measured by the first strain gauge 57 and the second strain gauge 59. The actually measured incident distorted wave is passed through a low-pass filter, and a high frequency of 10 kHz or more is removed. Furthermore, the time history of the incident distorted wave is subjected to zero correction in which the baseline value is zero. Each of the time axis strains in the first strain gauge 57 and the second strain gauge 59 obtained in this way is subjected to Fourier transform, and the frequency axis strain is obtained. A transfer function is derived from the frequency axis strain in the first strain gauge 57 and the second strain gauge 59. The ratio (X1: X2) between the distance X1 between the first strain gauge 57 and the rear end 53a of the input bar 53 and the distance X2 between the second strain gauge 59 and the rear end 53a of the input bar 53 is taken into consideration. Based on the transfer function, the frequency axis distortion at the rear end 53a of the input bar 53 is estimated. The frequency axis distortion is inversely Fourier transformed to obtain the time axis distortion (time history of distortion) εi of the incident distortion wave at the rear end 53a of the input bar 53.
[0077]
Similarly, the reflected strain wave reflected at the rear end 53 a of the input bar 53 and traveling toward the front end 53 b is actually measured by the second strain gauge 59 and the first strain gauge 57. From this actually measured reflected distortion wave, the time-axis distortion (distortion time history) εr of the reflected distortion wave at the rear end 53a of the input bar 53 is obtained.
[0078]
Further, the transmitted strain wave propagated to the output rod 55 through the test piece 70 is measured by the third strain gauge 61 and the fourth strain gauge 63 of the output rod 55, and the output rod 55 of the output rod 55 is measured from the measured transmitted strain wave. A time-axis strain (strain time history) εt of the transmitted strain wave at the front end 55a is obtained.
[0079]
From the obtained εi, εr, and εt, the strain rate ε ′ of the test piece 70 is calculated by the following formula 9.
(Formula 9)
ε ′ = (C0/ L) · (εi-εr-εt)
= ((E / ρ)1/2/ L) · (εi-εr-εt) --- (9)
(In Equation 9, C0Represents the propagation velocity (m / s) of the strain wave in the input bar and the output bar (stress bar), L represents the length (m) of the test piece, and E represents the longitudinal elastic modulus (N / m) of the stress bar.2) Is the density of the stress bar (kg / m)3Represents)
[0080]
Further, the strain ε of the test piece 70 is calculated from εi, εr, and εt by the following formula 10.
[0081]
(Formula 10)
Figure 0003910825
(In Equation 10, C0Represents the propagation velocity (m / s) of the strain wave in the input bar and the output bar (stress bar), L represents the length (m) of the test piece, and E represents the longitudinal elastic modulus (N / m) of the stress bar.2) Is the density of the stress bar (kg / m)3Represents)
[0082]
Further, the stress σ of the test piece 70 is calculated from εi, εr, and εt by the following formula 11.
(Formula 11)
σ = (E · A / (2As)) · (εi + εr + εt)
= (ED2/ (2 (Ds)2)) · (Εi + εr + εt) --- (11)
(In Formula 11, E is the longitudinal elastic modulus (N / m of a stress bar composed of an input bar and an output bar)2), And A is the cross-sectional area (m2As is the cross-sectional area of the test piece (m2), D represents the diameter (m) of the stress bar, and Ds represents the diameter (m) of the test piece)
[0083]
The strain time history of the test piece 70 thus obtained is shown in the graph of FIG. As shown in FIG. 13, the curve is gentle for a while after the peak P, but then becomes uneven. A point S at a gentle stage after the peak P is selected, a tangent to the curve at this point S is drawn, a relaxation time λ is derived from the intersection of this tangent and the time axis, and the curve obtained by the following equation 12 is By setting the curve after S, the entire strain time history can be a gentle curve (shown by a dotted line in FIG. 13). Thereby, the influence of noise on the material properties finally obtained can be removed.
(Formula 12)
ε (t) = ε0・ E-T / λ      --- (12)
(In Equation 12, ε0Represents the strain at the contact)
[0084]
Similarly, the entire stress time history can be made into a gentle curve by the following Equation 13, thereby removing the influence of noise on the material properties finally obtained.
(Formula 13)
σ (t) = σ0・ E-T / λ      ---- (13)
(In Equation (13), σ0Represents the stress at the contact)
[0085]
As described above, the strain time history and the stress time history of the test piece 70 are corrected.
[0086]
By the above method, strain time history data, strain rate time history data, and stress time history data at high speed and large deformation are obtained with a split Hopkinson bar tester.
[0087]
Hereinafter, the Example and comparative example of this invention are explained in full detail.
First, a golf ball was manufactured using a material mainly composed of urethane rubber. A material mainly composed of urethane was compression-molded with a die having a diameter of 42.8 mm at 160 ° C. for 30 minutes to form a one-piece ball.
[0088]
Using the above-mentioned split Hopkinson bar tester, the impact speed of the striking bar is 4 patterns of 7 m / s, 14 m / s, 20 m / s, and 25 m / s. (Room temperature 23 ° C., relative humidity 50%).
The maximum strain and maximum strain rate during measurement are shown below.
Impact speed 7m / s (maximum strain 0.12, maximum strain speed 1378 / s)
Impact speed 14m / s (maximum strain 0.24, maximum strain speed 2703 / s)
Impact velocity 20m / s (maximum strain 0.35, maximum strain rate 3898 / s)
Impact speed 25m / s (maximum strain 0.43, maximum strain speed 4716 / s)
Table 4 below shows the values of the phase angle δ at each collision speed because the measurement is performed with four patterns of collision speeds.
[0089]
[Table 1]
Figure 0003910825
[0090]
(Example)
Viscous resistance using strain, strain rate, and stress time history data measured with a split Hopkinson bar tester using a material mainly composed of urethane rubber as a test piece, and a viscoelastic model similar to that of the first embodiment. The simulation was performed considering the above. In the simulation method according to the first embodiment, the relationship between the strain, strain rate, and longitudinal elastic modulus is further input into the product model, thereby taking into account the change in the longitudinal elastic modulus due to the difference in strain and strain rate. The software used was a general-purpose software Ls-dyna940 manufactured by Japan Research Institute.
Table 1 shows the phase angle δ predicted by the simulation of the example.
[0091]
(Comparative example)
The time history data of strain, strain rate, and stress measured with a split Hopkinson bar tester using a material mainly composed of urethane rubber similar to that in the above example as a test piece was used. The longitudinal elastic modulus was variable, and a simulation was performed using strain under load, strain rate, and viscous resistance of the load. Except for the points described above, the example was the same as the example.
The phase angle δ predicted by the simulation of the comparative example is shown in Table 1 above.
[0092]
According to the finite element method analysis, a golf ball made of the above-mentioned material composed mainly of urethane rubber was subjected to three 200 g (same weight as the head) aluminum hollow rod models with initial speeds of 35, 40 and 45 m / s. The performance of the golf ball and the deformation state of the material when the golf ball was caused to collide with each other were simulated, and the coefficient of restitution in the analysis of the golf ball was calculated by the method described above. Table 2 below shows the coefficient of restitution in the analysis of Examples and Comparative Examples of the golf balls made of the material mainly composed of the urethane rubber obtained here.
[0093]
[Table 2]
Figure 0003910825
[0094]
Moreover, the value of the coefficient of restitution of the actual golf ball was measured by the following method by an experiment using the actual golf ball molded from the material mainly composed of the urethane rubber.
Table 2 shows the difference (%) between the coefficient of restitution in the analysis of Examples and Comparative Examples and the coefficient of restitution in the experiment using the following real objects.
[0095]
(Measurement of coefficient of restitution by experiment using actual golf ball)
As a method for measuring the coefficient of restitution of a golf ball, an aluminum hollow rod of 200 g (same as the weight of the head) is used as a substitute for a golf club head on a golf ball made of the above material at a room temperature of 23 ° C. Collisions were made at three speeds of 35, 40, and 45 m / s. The speeds of the hollow rod and the golf ball before and after the collision were measured, and the coefficient of restitution of the golf ball was calculated from the respective speed and weight.
The collision surface of the aluminum hollow rod was flat, and the front collision with the golf ball was flat. Since there is no angle on the collision surface like a golf club head, the golf ball does not rotate at the time of collision, so only the rebound performance of the golf ball can be evaluated.
[0096]
As shown in Table 1, the phase angle δ in the analysis of the example was almost the same as the phase angle δ of the experimental result at each collision speed. On the other hand, since the phase angle δ in the analysis of the comparative example is simulated using a load-to-load curve, the energy loss in the stress-strain curve is smaller than that in the example, and the comparative example at each collision speed. The value of the phase angle δ was smaller than that of the example, the value was smaller than the experimental result, and the difference was large. Therefore, it was confirmed that the example was able to simulate actual experimental results with high accuracy.
[0097]
In addition, as shown in Table 2, the value of the coefficient of restitution in the analysis of the examples is the value of the coefficient of restitution in an experiment actually obtained by actually producing a golf ball and any hollow rod speed. Also agreed well. The difference between the example and the actual experimental result was −4.82% to −5.98%, and it was confirmed that the actual experimental result was accurately simulated. On the other hand, the value of the coefficient of restitution in the analysis of the comparative example was small with respect to the experimental value because the value of the phase angle δ was small, and the difference was also large compared to the example. Specifically, the difference between the comparative example and the actual experimental result was 16.1% to 22.2%, which was greatly different from the experimental result.
[0098]
Based on the above, the difference in viscous resistance between the strain load state and the strain unloading or restoring state even when the strain is the same, depending on the case of the strain loading state and the strain unloading or restoring state. By performing the simulation of the present invention in consideration of the above, it was confirmed that the performance of the product made of the viscoelastic material whose material physical property is nonlinear can be predicted with high accuracy by the simulation.
[0099]
【The invention's effect】
As is clear from the above explanation, according to the present invention, in order to consider viscoelastic characteristics in the deviation component of strain and strain rate generated for each element, the strain and strain rate are decomposed into a volume component and a deviation component. Using the principal strain and principal strain rate of the deviation component converted to the principal strain coordinate system and principal strain rate coordinate system, the viscous resistance is determined for each axis of the principal strain coordinate system and principal strain rate coordinate system. The three components in the axial direction are variable for each element. For this reason, the analysis which considered the viscoelasticity of the viscoelastic material correctly can be performed about a three-dimensional direction, and the precision of a simulation can be improved.
[0100]
In addition, using a viscoelastic model that takes into account the viscosity of the viscoelastic material, the viscous resistance is derived separately for the strain loading state and the strain unloading or restoring state, and the strain, strain rate, and strain loading are derived. Because the simulation is performed by inputting the relationship of viscous resistance, which is different between the unloading state and the restoration state, into the product model, the actual material state can be accurately captured, and the nonlinearity of the viscoelastic material It is possible to accurately represent a phenomenon in which a physical property value indicating non-linearity changes nonlinearly depending on the deformation speed and size of a material.
[0101]
In addition, the strain, strain rate, and stress values used in the simulation are measured under conditions that assume the actual usage state of products made of viscoelastic materials, so they correspond to various deformation states of actual viscoelastic materials. Simulation can be performed.
[0102]
Therefore, even for products made of viscoelastic materials where the relationship between strain and strain rate changes depending on the deformation state of the material, and the material properties such as loss factor show nonlinearity, the performance and dynamic behavior of the product are determined by simulation. It is possible to analyze with high accuracy.
[0103]
Actually, the actual ball hit test can be simulated under the same strain and strain rate as when the golf ball was hit with a golf club head. Performance and deformation behavior can be predicted with high accuracy. This makes it easy to understand the material properties that determine the superiority or inferiority of the golf ball, which not only helps to improve the performance of the product, but also reduces the number of actual ball trials at the design stage of the golf ball. Time can be saved.
[Brief description of the drawings]
FIG. 1 is a diagram showing time history data of strain ε.
FIG. 2 is a diagram showing time history data of a strain rate ε ′.
FIG. 3 is a diagram showing time history data of stress σ.
FIG. 4 is a diagram showing a strain-stress diagram and a method of calculating a longitudinal elastic modulus.
FIG. 5 is a diagram showing a two-element Voigt model used as a viscoelastic model of the present embodiment.
FIG. 6 is a diagram showing time history data of viscous resistance.
FIG. 7 is a diagram showing the relationship between strain and viscous resistance.
FIG. 8 is a diagram showing the relationship between strain and strain rate.
FIG. 9 is a diagram showing a state of division by a mesh of a golf ball model.
FIGS. 10A and 10B show a collision situation of an aluminum hollow rod model with a golf ball model. FIG. 10A is a diagram before the collision, FIG. 10B is a diagram during the collision, and FIG.
FIG. 11 is a diagram illustrating a method for calculating a loss factor.
FIG. 12 is a schematic front view showing a split Hopkinson bar testing machine.
FIG. 13 is a diagram showing a state of strain time history of a test piece.
[Explanation of symbols]
10 Golf club model
11 elements
12 nodes
20 Aluminum hollow bar model

Claims (6)

粘弾性材料からなる製品の実使用状態を想定した測定条件下で、上記粘弾性材料に生じるひずみ、ひずみ速度、応力の時々刻々の値を測定し、
上記ひずみ、ひずみ速度、応力の時刻歴データと、上記粘弾性材料の粘性を考慮した粘弾性モデルとから、ひずみの負荷状態時とひずみの除荷あるいは復元状態時とで場合分けして上記粘弾性材料の粘性抵抗の時刻歴データを導出し、
上記粘弾性材料からなる製品を、解析対象の製品モデルとして設定し、該製品モデルを多数の要素に分割し、上記ひずみ、ひずみ速度、粘性抵抗の関係を入力し、上記ひずみの負荷状態時とひずみの除荷あるいは復元状態時とでの粘性抵抗の違いを考慮して有限要素法により解析を行い、上記粘弾性材料からなる製品モデルの性能を予測する際に、
上記要素毎に生じる全体座標系のひずみ及びひずみ速度を偏差成分と体積成分に分解し、上記偏差成分のひずみ及びひずみ速度を、全体座標系から主ひずみ座標系及び主ひずみ速度座標系に変換し、
上記変換した偏差主ひずみ及び偏差主ひずみ速度を用い、主ひずみ座標系及び主ひずみ速度座標系の各座標軸毎に、ひずみが同じ値でもひずみの負荷状態時とひずみの除荷あるいは復元状態時とで異なった粘性抵抗を決定していることを特徴とする粘弾性材料からなる製品の性能予測のためのシミュレーション方法。
Under the measurement conditions that assume the actual use state of the product made of viscoelastic material, measure the momentary value of strain, strain rate, stress generated in the viscoelastic material,
Based on the time history data of the strain, strain rate, and stress, and the viscoelastic model taking into account the viscosity of the viscoelastic material, the viscoelastic material is classified into cases when the strain is loaded and when the strain is unloaded or restored. Deriving time history data of viscous resistance of elastic materials,
The product made of the viscoelastic material is set as a product model to be analyzed, the product model is divided into a number of elements, the relationship between the strain, strain rate, and viscous resistance is input, and when the strain is loaded When analyzing the finite element method in consideration of the difference in viscous resistance between unloading or restoring the strain and predicting the performance of the product model made of the above viscoelastic material,
The strain and strain rate of the global coordinate system generated for each element are decomposed into a deviation component and a volume component, and the strain and strain rate of the deviation component are converted from the global coordinate system to the main strain coordinate system and the main strain rate coordinate system. ,
For each coordinate axis of the main strain coordinate system and the main strain rate coordinate system, using the converted deviation main strain and deviation main strain rate, even when the strain is the same, the strain is loaded and the strain is unloaded or restored. A simulation method for predicting the performance of a product made of a viscoelastic material, characterized in that different viscous resistances are determined in
上記要素毎に主ひずみの大きさであるノルムを計算し、シミュレーション時の前ステップのノルムと比較することによって、ノルムが増加している場合は負荷、ノルムが減少している場合は除荷と判断し、ひずみの負荷と除荷の状態を判定している請求項1に記載の粘弾性材料からなる製品の性能予測のためのシミュレーション方法。Calculate the norm, which is the magnitude of the principal strain, for each of the above elements and compare it with the norm of the previous step at the time of simulation, so that if the norm is increasing, it is loaded, and if the norm is decreasing, unloading is performed. The simulation method for predicting the performance of a product made of a viscoelastic material according to claim 1, wherein the determination is made and the state of strain loading and unloading is determined. 上記粘弾性材料は、材料物性が非線形性を示す材料である請求項1または請求項2に記載の粘弾性材料からなる製品の性能予測のためのシミュレーション方法。The simulation method for predicting the performance of a product made of a viscoelastic material according to claim 1, wherein the viscoelastic material is a material whose material physical properties are nonlinear. 上記ひずみ、ひずみ速度、応力をスプリットホプキンソン棒試験機により測定している請求項1乃至請求項3のいずれか1項に記載の粘弾性材料からなる製品の性能予測のためのシミュレーション方法。The simulation method for predicting the performance of a product made of a viscoelastic material according to any one of claims 1 to 3, wherein the strain, strain rate, and stress are measured by a split Hopkinson bar tester. 上記測定時に粘弾性材料に生じる、上記ひずみの最大値が0.05〜0.50の範囲であり、あるいは/及び、上記ひずみ速度の最大値が500/s〜10000/sの範囲である請求項1乃至請求項4のいずれか1項に記載の粘弾性材料からなる製品の性能予測のためのシミュレーション方法。The maximum value of the strain generated in the viscoelastic material during the measurement is in the range of 0.05 to 0.50, and / or the maximum value of the strain rate is in the range of 500 / s to 10000 / s. A simulation method for predicting the performance of a product made of the viscoelastic material according to any one of claims 1 to 4. 上記粘弾性材料はゴルフボール用材料とし、上記製品モデルをゴルフボールとし、
上記ゴルフボールと、ゴルフクラブヘッドを想定した打撃物との衝突現象をシミュレーションし、該衝突時のゴルフボールの挙動を予測している請求項1乃至請求項5のいずれか1項に記載の粘弾性材料からなる製品の性能予測のためのシミュレーション方法。
The viscoelastic material is a golf ball material, the product model is a golf ball,
The viscosity according to any one of claims 1 to 5, wherein a collision phenomenon between the golf ball and a hitting object assuming a golf club head is simulated to predict the behavior of the golf ball at the time of the collision. A simulation method for predicting the performance of products made of elastic materials.
JP2001341117A 2001-07-18 2001-11-06 Simulation method for performance prediction of products made of viscoelastic materials Expired - Fee Related JP3910825B2 (en)

Priority Applications (2)

Application Number Priority Date Filing Date Title
JP2001341117A JP3910825B2 (en) 2001-11-06 2001-11-06 Simulation method for performance prediction of products made of viscoelastic materials
US10/481,754 US7130748B2 (en) 2001-07-18 2002-07-16 Simulation method for estimating performance of product made of viscoelastic material

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
JP2001341117A JP3910825B2 (en) 2001-11-06 2001-11-06 Simulation method for performance prediction of products made of viscoelastic materials

Publications (2)

Publication Number Publication Date
JP2003139668A JP2003139668A (en) 2003-05-14
JP3910825B2 true JP3910825B2 (en) 2007-04-25

Family

ID=19155205

Family Applications (1)

Application Number Title Priority Date Filing Date
JP2001341117A Expired - Fee Related JP3910825B2 (en) 2001-07-18 2001-11-06 Simulation method for performance prediction of products made of viscoelastic materials

Country Status (1)

Country Link
JP (1) JP3910825B2 (en)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US6829563B2 (en) * 2001-05-31 2004-12-07 Sumitomo Rubber Industries, Ltd. Simulation method for estimating performance of product made of viscoelastic material
JP5039178B2 (en) * 2010-06-21 2012-10-03 住友ゴム工業株式会社 Method for simulating products made of viscoelastic materials
CN113109189B (en) * 2021-03-26 2023-07-11 北京工业大学 Method for determining cyclic stress strain of resin matrix composite material by considering frequency

Also Published As

Publication number Publication date
JP2003139668A (en) 2003-05-14

Similar Documents

Publication Publication Date Title
US7254492B2 (en) Method of computing energy loss generated in viscoelastic material and method for evaluating energy loss of golf ball by using method of computing energy loss
US6829563B2 (en) Simulation method for estimating performance of product made of viscoelastic material
EP2642417B1 (en) Method for estimating energy loss of high polymer material
US7130748B2 (en) Simulation method for estimating performance of product made of viscoelastic material
JP3466584B2 (en) Simulation method for predicting the performance of products made of viscoelastic materials
Karagiozova et al. Impact of aircraft rubber tyre fragments on aluminium alloy plates: II—Numerical simulation using LS-DYNA
JP3910825B2 (en) Simulation method for performance prediction of products made of viscoelastic materials
JP3910826B2 (en) Simulation method for performance prediction of products made of viscoelastic materials
JP3466585B2 (en) Simulation method for predicting the performance of products made of viscoelastic materials
JP3466590B2 (en) Simulation method for predicting the performance of products made of viscoelastic materials
JP4223178B2 (en) Ball rebound characteristics prediction method
JP3466583B2 (en) Simulation method for predicting the performance of products made of viscoelastic materials
Daiyan et al. Numerical simulation of low-velocity impact loading of a ductile polymer material
EP2397956B1 (en) Method for simulating article made of viscoelastic material
US6671642B2 (en) Method of evaluating energy loss of golf ball
JP3597517B2 (en) Golf ball energy loss evaluation method
Tamaogi et al. Identification of mechanical models for golf ball materials using a viscoelastic split hopkinson pressure bar
Gan et al. Effect of Rubber Strain Rate on Performance of Aircraft Tire
Mase T. Mase and AM Kersten,” Experimental evaluation of a 3-D hyperelastic, rate-dependent golf ball constitutive model,” The Engineering of Sport 5-Conference Proceedings of, Vol. 2, pp. 238-244, 2004
KR101560060B1 (en) Appratus for testing profile of a tire and tire performance prediction system
JP2006275926A (en) Simulation method capable of predicting performance of product
Afdhal et al. Numerical simulation of SHPB to measure the mechanical properties of aluminium foam material at high strain rate by using MAT 163 modified crushable foam
Makowski Determination of Material Properties of Basketball. Finite Element Modeling of the Basketball Bounce Height
Han et al. Development of a Flex-PLI System Model and Investigations of Injury
Liu et al. Golf ball impact: material characterization and transient simulation

Legal Events

Date Code Title Description
A621 Written request for application examination

Free format text: JAPANESE INTERMEDIATE CODE: A621

Effective date: 20040916

A711 Notification of change in applicant

Free format text: JAPANESE INTERMEDIATE CODE: A711

Effective date: 20050520

RD02 Notification of acceptance of power of attorney

Free format text: JAPANESE INTERMEDIATE CODE: A7422

Effective date: 20050608

TRDD Decision of grant or rejection written
A01 Written decision to grant a patent or to grant a registration (utility model)

Free format text: JAPANESE INTERMEDIATE CODE: A01

Effective date: 20070123

A61 First payment of annual fees (during grant procedure)

Free format text: JAPANESE INTERMEDIATE CODE: A61

Effective date: 20070125

R150 Certificate of patent or registration of utility model

Free format text: JAPANESE INTERMEDIATE CODE: R150

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20100202

Year of fee payment: 3

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20110202

Year of fee payment: 4

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20120202

Year of fee payment: 5

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20120202

Year of fee payment: 5

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20130202

Year of fee payment: 6

FPAY Renewal fee payment (event date is renewal date of database)

Free format text: PAYMENT UNTIL: 20140202

Year of fee payment: 7

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

R250 Receipt of annual fees

Free format text: JAPANESE INTERMEDIATE CODE: R250

LAPS Cancellation because of no payment of annual fees