6.1 接触解析の基本(摩擦あり)

トップページは<a href="https://qiita.com/Sagittarius_Chiron/items/fd806361ceca4864bc79">こちら</a>です。
1.0、2.0、4.0の記事の読了を前提として、操作についての記載を省略している場合があります。
ver1.1.1では近接している面を自動認識する機能が追加されていますが、ここでは使用しません。

### 解析内容
<dl>
  <dt>使用ファイル</dt>
  <dd>6_01_Base.stp と 6_01_Top.stp　　(注1)</dd>
  <dt>単位</dt>
  <dd>mm-ton</dd>
  <dt>材料物性値:Top</dt>
  <dd>Aluminum　　(注2)</dd>
  <dd>縦弾性係数：E=70600 [MPa]</dd>
  <dd>ポアソン比：ν=0.343</dd>
  <dt>材料物性値:Base</dt>
  <dd>Steel　　(注2)</dd>
  <dd>縦弾性係数：E=212000 [MPa]</dd>
  <dd>ポアソン比：ν=0.293</dd>
</dl>
荷重状態
　load(面)：変位 Dy = - 0.5, Dx = Dz = 0 [mm]
　fix(面)　：変位 Dx =0.5, Dy = Dz = 0 [mm]
　symm(面)：対称面
　TopとBaseの間の初期隙間 = 0 [mm]
　接触の解法＝ペナルティ法
　摩擦係数　μ=1.5
　法線方向のバネ定数　E_N = 70000 [N/mm³]
　接線方向のバネ定数　E_T = 7000 [N/mm³]
<img width="75%" alt="Fig611.PNG" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888272/16412084-a0b4-84b2-6d9c-5ea585f9defb.png">

<ul>
<li>(注1)</li>
元の例題の形状セットには含まれていないため作成が必要でした。<br>
作成したファイルを<a href="https://drive.google.com/drive/u/2/folders/1j3PD-Ks2W_dTvEu-m0jNoZzfGDIQQz1e"><strong>こちら</strong></a>に置きました。<br>

<li>(注2)</li>
この値は元の例題での値です。ここでは PrePoMax の Material Library に登録済の近い値の材料を使います。
</ul>

### 形状の読込
PrePoMaxの新規起動または新規作業開始をしてください。
unit system:<b>mm,ton,s,℃</b> に設定してください。

6_01_Top.stp をインポートしてください。
オブジェクト・ブラウザのGeometryタブのツリーにSolid_part-1が作成され、グラフィック・パネルには部品が表示されます。

6_01_Base.stpをインポートしてください。
オブジェクト・ブラウザのGeometryタブのツリーにSolid_part-2が作成され、グラフィック・パネルには部品が表示されます。

ツリーに表示される名前の変更ができますが、ここでは
・Top部品のインデックス = 1
・Base部品のインデックス = 2
という規則で設定していきます。

### メッシュ作成
Solid_part-1にCreate Meshを実行してメッシュを生成してください。
Solid_part-2にCreate Meshを実行してメッシュを生成してください。
※　番号順に操作してください。

### 材料物性値の設定

Material Library から材料物性値を呼び出してください。

　FE ModelタブのツリーのMaterialsの右クリックメニューのMaterial Library。
　メニューバーの場合はProperty→Material library。

　Steel：S185
　Aluminum：1060_Alloy
<img width="50%" alt="Fig402.PNG" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888272/55433bcf-7a43-5db9-849d-5d4ec7fe88c4.png">
　ツリーのMaterialsの下に S185 と 1060_Alloy が表示されます。

セクションの設定をしてください。

　FE ModelタブのツリーのSectionsの右クリックメニューのCreate。
　メニューバーの場合はProperty→Section→Create。

　<b>Type</b> = Solid selection
　<b>Properties</b>
　Name = Solid_selection-1 (デフォルトのまま)
　Material =  1060_Alloy
　Region type = Part name
　Part = Solid_part-1

　「OK-New」を左クリック

　<b>Type</b> = Solid selection
　<b>Properties</b>
　Name = Solid_selection-2 (デフォルトのまま)
　Material =  S185
　Region type = Part name
　Part = Solid_part-2

　「OK」を左クリック
<img width="75%" alt="Fig403.PNG" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888272/3cb6a982-d4cf-509b-f473-b8dc4b0a4cfa.png">


ツリーのSectionsの下に Solid_selection-1 と Solid_selection-2 が表示されます。

### 面の命名
　面に名前を付けてください。

　FE ModelタブのツリーのMeshのSurfaces の右クリックメニューのCreate。
　メニューバーの場合はModel→Surfce→Create

1) <b>Properties</b>
　 Name = fix
　 Surface type = Element
　 Region type = Selection
2) Surfaces, edges and vertices
3) 面を選択してください。
4) 「Ok-New」または「OK」を左クリックしてください。

同様にして、あと４つの面の設定をしてください。

Name = contBase
Name = contTop
Name = load
Name = symm

contTopを選択する時は、下側のBase部品を非表示状態にしたほうがやりやすいでしょう。非表示と表示を切り替える方法は下記です。
<img width="70%" alt="Fig509.PNG" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888272/52b646ff-e2e9-b3ed-759b-ef481b034fc4.png">
symm(面)を設定するときは２つの面を選択する必要があります。複数の面を選択するときは、Shittキーを押しながら左クリックすると面を追加できます。CTRLキーを押しながら左クリックすると選択済みの面をセットから外すことができます。

作業が終わったあとは、ツリーのSurfacesの下に fix、contBase、contTop、load、symm が表示されているはずです。

### 接触面(相互作用)の設定
接触面の設定をしてください。

　FE ModelタブのツリーのMeshのContact → Surfaces interactions の右クリックメニューのCreate。
　メニューバーの場合はInteraction → Contact → Surfaces interaction → Create


1)  <b>Data</b>
　　Interaction name = Surface_interaction-1(デフォルトのまま)
2)　Interaction models
　　Selected = Surface behavior
　<b>Properties</b>
　　Data
　　Type = Linear
　　K = 10000000 N/mm³(デフォルト)から 70000 N/mm³に変更する
　　σ infinity = 1 MPa(デフォルトのまま)
　　Co = Default(デフォルトのまま)

3)　Interaction models
　　Selected = Friction
　<b>Properties</b>
　　Data
　　Friction coefficient = 1.5
　　Stick slope = Default
4) 「OK」を左クリックしてください。

<img width="100%" alt="Fig612.PNG" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888272/26069a71-3abc-9718-abb7-5f3a7bb3c8f2.png">

≪補足≫
CaliculiXのマニュアルを見ると、以下のことが書いてありました。
「Type = Hard」という設定は「Type = Linear」のデフォルト条件である。
「Type = Linear」のKは、法線方向の接触剛性である。
「Type = Linear」のSlick slope(λ)と、接線方向の接触剛性Ktの関係式は下記である。
　　Kt = λ × Ap
　　Apは“スレーブ面のノードの代表”らしいのですが、意味がよくわかりません。

つまり、法線方向の接触剛性Kは、元の例題の E_N のことです。
元の例題の E_T は PrePoMaxのGUI画面では設定できません。Slick slope(λ)の指定によってE_Tを指定することになるのですがApの意味が判らないため、元の例題のE_T=7000に合わせた換算はしないでデフォルトのままとしておきます。CaliculiXの*.inpファイルをエクスポートすると、Slick slope(λ)は未指定の空欄(つまりデフォルト値)となっていました。したがって、この設定で元の例題に相当しているとみなします。

接触面ペアの設定をしてください。
　FE ModelタブのツリーのMeshのContact → Contact pairs の右クリックメニューのCreate。
　メニューバーの場合はInteraction → Contact → Contact pair → Create

1) <b>Data</b>
　　Name = Contact_Pair-1 (デフォルトのまま)
　　surface interaction = Surface_interaction-1 (デフォルトのまま)
　　Method = Surface to surface (デフォルトのまま)
　　Adjust = Yes (デフォルトのまま)
　　Adjustment size = Default (デフォルトのまま)
　<b>Master Region</b>
　　Master region type = Surface name
　　Master surface = contBase
　<b>Slave Region</b>
　　Slave region type = Surface name
　　Slave surface = contTop
3) 「Ok」を左クリックしてください。
<img width="50%" alt="Fig613.PNG" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888272/4bb9bb14-4bac-4465-1ae2-4687954751e5.png">

### 解析設定
解析ステップを新規作成してください。

　FE ModelタブのツリーのStepsの右クリックメニューのCreate。
　メニューバーの場合はStep→Step→Createです。

　<b>Type</b> = Static step
　<b>Data</b>
　Name = Step-1 (デフォルトのまま)
　Solver = Default (デフォルトのまま)
　Nlgeom = <font color="Red">On</font> 
　Incrementation = Default (デフォルトのまま)
<img width="200" alt="Fig08.PNG" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888272/28f439ba-53a5-b205-174c-dd051bbd162a.png">

固定する面を作成してください。

　FE ModelタブのツリーのStep-1/BCsの右クリックメニューのCreate。
　メニューバーの場合はStep→BC→Createです。

　<b>Type</b> = Displacement/Rotation
　<b>Properties</b>
　Name = Displacement_rotation-1 (デフォルトのまま)
　Region type = Surface name
　Surfaces = symm
　U1 = Unconstrained
　U2 = Unconstrained
　U3 = 0 mm
　UR1 = Unconstrained
　UR2 = Unconstrained
　UR3 = Unconstrained
　Color = Lime (デフォルトのまま)
　<br>
　<b>Type</b> = Displacement/Rotation
　<b>Properties</b>
　Name = Displacement_rotation-2 (デフォルトのまま)
　Region type = Surface name
　Surfaces = load
　U1 = 0 mm
　U2 = -0.5 mm
　U3 = 0 mm
　UR1 = Unconstrained
　UR2 = Unconstrained
　UR3 = Unconstrained
　Color = Lime (デフォルトのまま)

  <br>
　<b>Type</b> = Displacement/Rotation
　<b>Properties</b>
　Name = Displacement_rotation-3 (デフォルトのまま)
　Region type = Surface name
　Surfaces = fix
　U1 = 0.5 mm
　U2 = 0 mm
　U3 = 0 mm
　UR1 = Unconstrained
　UR2 = Unconstrained
　UR3 = Unconstrained
　Color = Lime (デフォルトのまま)

≪補足≫
　Displacement_rotation-2 がTop部品を押さえつける条件で、
　Displacement_rotation-3 がBase部品を水平移動させる条件です。

　今回は、静的計算(Static step)で、Nlgeom(幾何非線形)=On という設定です。

　この設定で計算すると、２つの条件が同時進行で
　Top部品の移動量　(　0,　0,　0　)　→　(　0,　-0.5,　0　)
　Base部品の移動量　(　0,　0,　0　)　→　(　0.5,　0,　0　)
　という変化をする計算になります。

　元の例題は、Displacement_rotation-2が最初から (　0,　-0.5,　0　) という状態で、Displacement_rotation-3がランプ関数にしたがい(　0,　0,　0　)　→　(　0.5,　0,　0　)と変化する計算です。
　元の例題に倣った設定での計算は、あとで行います。

　これらの３つの違いを示した模式図を下に示します。
<img width="80%" alt="Fig614.PNG" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888272/def7210c-aa2e-fdc6-64ad-8669ad7abdbd.png">


≪ファイル保存≫
いったんファイルを保存してください。名前は任意ですが、ここでは06_10フォルダーの06_10_a.pmxファイルとしました。

### ソルバー計算
ソルバー計算を実施してください。

表示例を示します。変形形状は True scale です。
<img width="75%" alt="Fig615.PNG" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888272/4c508620-e61f-6a8a-840c-04342a361ff7.png">
<img width="75%" alt="Fig616.PNG" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888272/ad8823db-e704-abb4-2805-19e79ea02be6.png">


≪補足≫
参考にさせていただいた元の例題(SalomeMeca)を熟読した人から質問があるかもしれないので、補足します。
元の例題には、『接触問題では２次要素は収束し難く計算が遅いため線形要素を使ったほうが良い』という記載があります。CalculiXのマニュアルには、『線形要素を使うときは接触面の計算方法(Method)を Node to surfaceで良いが、２次要素を使うときはMethodを Surface to surfaceとする。』旨の記載があります。今回は「Method = Surface to surface」で計算しています。

≪補足≫
元の例題の(2)では、ペナルティ法に代えてラグランジェ法（ラグランジュの未定乗数法）での計算をしています。CalculiXのマニュアルを読むとペナルティ法しかできないようです。自分はCaliculiXの初心者ですので、もしも勘違いがあれば御指摘ください。

### ２段階計算

ファイルを上書き保存してください。
そして、別名で保存してください。ファイル名は06_10フォルダーの06_10_b.pmxとしました。

解析ステップを追加作成してください。

　FE ModelタブのツリーのStepsの右クリックメニューのCreate。
　メニューバーの場合はStep→Step→Createです。

　<b>Type</b> = Static step
　<b>Data</b>
　Name = Step-2 (デフォルトのまま)
　Solver = Default (デフォルトのまま)
　Nlgeom = <font color="Red">On</font> 
　Incrementation = Default (デフォルトのまま)

下の図のように、Step-1の次にStep-2が追加されます。Step-1の境界条件はStep-2に引き継がれます。
<img width="50%" alt="Fig617.PNG" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888272/1d743aae-aa4d-97e2-c5c9-b87f306e5c98.png">

<b>Step-1 </b>の<b> Displacement_rotation-3 </b>
　Surfaces = fix
　U1 = 0 mm  に変える
　U2 = 0 mm
　U3 = 0 mm
　UR1 = Unconstrained
　UR2 = Unconstrained
　UR3 = Unconstrained

ソルバー計算を実施してください。

表示例を示します。変形形状は True scale です。さきほどの結果と同じです。
<img width="75%" alt="Fig618.PNG" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888272/a0f60cd0-369f-b612-8817-6a944993fbd4.png">
<img width="75%" alt="Fig619.PNG" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888272/f0bf71dc-3ba1-7256-4195-71349e61c393.png">

##### ２段階計算のバグについて
この記事を書いたあとに、２段階計算で一部の条件が引き継がれないというPrePoMaxの<b>“仕様(笑)”</b>が判りました。
発見に至った経緯は「5.0 線形熱応力」の最後に記載してあります。Initial Condition で設定した条件がソルバー計算の最初で初期化され、２段階計算のときは前のステップの境界条件をリセットされます。特に、前のステップの境界条件がリセットされるというは、非定常計算をする時には完全にダメです。
この記事では、２段階計算をしています。
本来ならば、第１ステップでTop部品へ印加したのちに第２ステップはその状態を保持します。“仕様(笑)”のせいでリセットされてしまいます。したがって、下の図の「この例題」と同じ状態です。
<img width="80%" alt="Fig614.PNG" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/888272/def7210c-aa2e-fdc6-64ad-8669ad7abdbd.png">
</br>
(他の記事は、2022/02/04 22:55頃に修正済です。この記事の追記修正を忘れていました。)