円筒のはめあい

今週は、アセンブリ構造解析の事例として、部品の接触解析を行いました。  
解析は、焼き嵌めを想定し2円筒がはめあったときの部品の強度を計算しました。  
2022.1.8更新　機械便覧の計算式を用いた円周応力計算をpythonで行った結果を追加しました。  

使用ソフト：SalomeMeca2019+Ubuntu18.04(with VirtualBox)  
　　　　　　SalomeMeca2019+CAELinux2020Lite(with VirtualBox)

==2023.2.26 SalomeMeca2021 for windows(windows11)==でも動作確認しました。

## １．解析概要

２つの円筒がはめあったときにおける外側円筒（円柱２）の円周応力を計算します。  



（１）モデル  
上図のように、軸対称モデルにしており、はめあい面に接触が発生すっることから非線形解析による接触解析を行いました。  
・円柱1：内径10mm x 外径20mm x 高さ10ｍｍ  
・円柱2：内径19.98mm x 外径30mm x 高さ11ｍｍ  
・しめしろ： 0.02mm  
・ヤング率：200000MPa、ポアソン比0.3（鉄系材料）  
・出力：外輪はめあい面の円周応力  
・メッシュ：四角形１次要素  
・接触に関連するパラメータ：SalomeMecaの初期値

（2）要素  
　4面体1次要素（自動メッシュ）  
　メッシュサイズは、2条件設定しました。  
　＜条件１＞  
　径方向（X方向）：5分割、高さ方向（Y方向）：10分割  
　＜条件２＞  
　径方向（X方向）：20分割、高さ方向（Y方向）：100分割

## ２．解析結果

結果を以下に示します。  
円周応力はZ方向であり、紙面に垂直な方向です。

![picture.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782081/13fb511a-e7ac-f450-286d-f4b8053043de.png)


##   
３．検算

## ３．１　Webツールによる公式による計算

はめあい時の応力は、材料力学の公式で算出可能です。  
計算が少々複雑なので、インターネットにあるWeb計算ツールを活用しまして、算出すると121.9MPaでした。  
SalomeMecaの応力値は、条件１は124.8MPa、条件２は122.2であることから、材料力学の公式による計算値とほぼ一致していました。  

> Web計算ツール：HP「計算技術研究所」内の圧入計算より
> 
> URL:[https://gijyutsu-keisan.com/webapp/mech/calc_pressfit/explain.php](https://gijyutsu-keisan.com/webapp/mech/calc_pressfit/explain.php)

## ３．２　機械便覧による公式による計算

材料力学等で円筒の円周応力による計算式があります。今回は、機械便覧による公式を用いて作成を行いました。

>   
> はめあい面の面圧Pmの計算式
> 
> $$Pm=\frac{\frac{\delta}{c}+2( \frac{Pa}{E1}* \frac {a^2}{c^2-a^2}+ \frac{Pb}{E2}* \frac {b^2}{b^2-c^2} )}{ (\frac{1-\nu1}{E1}+ \frac{1-\nu2}{E2})+2( \frac{1}{E1}* \frac {a^2}{c^2-a^2}+ \frac{1}{E2}* \frac {b^2}{b^2-c^2} ) }$$
> 
> Pm:はめあい面の面圧、σθ:半径方向位置rにおける円周応力  
> a:内輪の内径半径、b:外輪の外径半径、c:内輪の外径内径、δ:しめしろ（半径分）、E1:内輪のヤング率、E2:外輪のヤング率、v1:内輪のポアソン比、v2:外輪のポアソン比、Pa：内輪内径にかかる面圧、Pb：円筒外輪外径にかかる面圧  
> R:半径寸法、k=b/a、R=r/a  
> 
> 引用文献：機械工学便覧 α基礎編　第6章　圧力・遠心力などを受ける円板、円筒、玉より  
> 　　　　　σθの計算式：α3-63（表6.1　内輪及び外圧を受ける円筒の弾性応力より）  
> 商品ページ(amazon)：[こちら](https://amzn.to/3t8KJlO)

>   
> 円筒の内径に面圧がかかったときの任意の半径位置rにおける円周応力 σθの計算式
> 
> $$\sigma_{\theta} = -\frac{\frac{k^2}{R^2+1}}{k^2-1}Pa$$
> 
> σθ:半径方向位置rにおける円周応力  
> a:円筒の内径半径、b:円筒の外径半径、Pa：円筒内径にかかる面圧  
> r:半径寸法、k=b/a、R=r/a  
> 
> 引用文献：機械工学便覧 α基礎編　第6章　圧力・遠心力などを受ける円板、円筒、玉より  
> 　　　　　Pmの計算式： α3-65、  
> 　　　　　σθの計算式：α3-63（表6.1　内輪及び外圧を受ける円筒の弾性応力より）  
> 商品ページ(amazon)：[こちら](https://amzn.to/3t8KJlO)

これを計算するのは間違えそうでしたので、google Colabを用いて計算をしました。  
書式はjupiter notebookと同じため、jupiternotebookでも同様の計算ができます。  
pythonをにはsympyという代数計算ライブラリがあり、計算項が多い計算式をつかええます。

1)ライブラリの宣言

```
from sympy import *
from sympy.abc import *
from sympy import I, pi, E
import numpy as np
import matplotlib.pyplot as plt
```

  
2)変数の宣言

```
a =Symbol('a')#内輪の内径半径
b= Symbol('b')#外輪の外径半径
c= Symbol('c')#はめあい面
r= Symbol('r')#半径方向座標
delta= Symbol('delta')#しめしろ
E1= Symbol('E1')#内輪ヤング率
E2= Symbol('E2')#外輪ヤング率
v1= Symbol('v1')#内輪ポアソン比
v2= Symbol('v2')#外輪ポアソン比
Pa= Symbol('Pa')#内輪にかかる内圧
Pb= Symbol('Pb')#外輪にかかる外圧
E= Symbol('E')#ヤング率
v= Symbol('v')#ポアソン比
epsiron_r = Symbol('epsiron_r') #半径方向ひずみ
epsiron_siita = Symbol('epsiron_siita')#円周方向ひずみ
sigma_r =Symbol('sigma_r')#面圧
sigma_siita =Symbol('sigma_siita')#円周応力
pm0=Symbol('pm0')#外輪にかかる外圧
```

3)はめあい面の面圧pmの計算式  
今回は外圧Pb=0、内圧Pa=0のため、分子の式： pm1_numeratorはdelta/cのみです。  
また、内輪、外輪ともに同じ材料を使っているため、E=E1=E2、v=v1=v2としています。  
計算すると以下の計算式が表示されます。

```
pm2_numerator = delta/c
pm2_denomirator =(1-v)/E+(1+v)/E+2*(1/E*a**2/(c**2-a**2)+1/E*c**2/(b**2-c**2))
print(pm_numerator)
print(pm_denomirator)
pm2 = pm2_numerator/pm2_denomirator
pm2
```

![](https://aytechlab.com/wp-content/uploads/2022/01/image-3.png)

4)諸元の代入による面圧の計算

```
#面圧の計算
pm2value=pm2.subs([(a,5),(b,15),(c,10),(delta,0.01),(Pa,0),(Pb,0),(E,200000),(v,0.3)])
pm2value
```

5)はめあい面の円周応力の 計算式  
外輪の場合、 r=b/c、r=cと計算しました。このとき外径面は 外輪内径と内輪外径の差は0.05％しか差がないため、外輪内径を内輪外径として近似しています。

```
sigma_siita =(b**2/c**2/(r**2/c**2)+1)/(b**2/c**2-1)*pm0
sigma_siita
```

![](https://aytechlab.com/wp-content/uploads/2022/01/image-5.png)

6) 諸元の代入によるはめあい面の円周応力の計算

```
#円周応力
sigma_siita_value =sigma_siita.subs([(b,15),(c,10),(r,10),(pm0,pm2value)])
sigma_siita_value
```

計算式は、121.88MPaとなり、Web計算ツールと同じとなります。

7)外輪の円周応力の分布  
matplotlibで外輪の半径方向と円周応力をプロットさせると以下のグラフができます。  
outeringは、外輪内径から外輪外径までの半径分、stressは円周応力です。

![image-4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3782081/ad12a43c-6716-e088-92ed-7b3ab799218e.png)


## ４．まとめ

SalomeMecaのはめあい計算は、材料力学の公式とほぼ一致していました。  
そのため、断面形状が不均一な場合な円筒でも、円筒にかかる応力計算を推定できると考えます。  
  
なお、この計算は、2020.12.13に開催された第20回オープンCAE勉強会@関東（構造など）でも、発表したところ助言をもらいました。  
・接触解析は、パラメータが複雑で収束がしにくい問題で難しい。  
・比較的簡単な接触解析のため、SalomeMecaの初期設定値でもエラーがでないで解析ができた。  
また、簡易的な検算ツールを公開いただきました技術計算製作所様にはこの場にてお礼申し上げます。  
  
pythonと google Colab またはjupiternotebookを用いた検算は、数式を見ながら計算ができるため、本記事で計算したはめあいのような計算項が多い計算を数式を確認しながら計算できるため、ミスなく計算ができます。  

＜参考資料＞  
*.comm：解析ファイル、*.med：メッシュファイル  
メッシュファイルは、条件１のみ添付します。

[Pressfit.comm](https://aytechlab.com/wp-content/uploads/2021/03/Pressfit.comm)[ダウンロード](https://aytechlab.com/wp-content/uploads/2021/03/Pressfit.comm)

[Compound_Mesh_1.med](https://aytechlab.com/wp-content/uploads/2021/03/Compound_Mesh_1.med)[ダウンロード](https://aytechlab.com/wp-content/uploads/2021/03/Compound_Mesh_1.med)

## ５．書籍化しました！

この解析の詳細は、技術書典12にて書籍の形で作成しました。詳細は本をお読みいただければと思います。  
書籍名：自宅で始める応力解析２ -SalomeMeca2019で行う始める円筒のはめあい解析-  
URL:[https://aytechlab.com/pr/202201bookfem/](https://aytechlab.com/pr/202201bookfem/)