Loading
接触面平滑化機能（Nagataパッチ）および接触反力ダンピング機能の追加
Close #747 (closed) Close #748 (closed)
概要
曲面転がり接触における数値的振動・収束性悪化に対処するため、以下の2機能を追加する。
- 接触面平滑化機能（Nagataパッチ補間） —
!CONTACTにSMOOTHING=NAGATAオプションを追加 - 接触反力ダンピング機能 —
!CONTACTにDAMP_ALPHA・DAMP_GACTパラメータを追加
1. 接触面平滑化機能（Nagataパッチ）
問題設定
有限要素の接触面は多角形要素の平面集合であるため、曲面接触では要素面法線の不連続が接触力の振動を招く。Nagataパッチ法は各要素の頂点法線を用いて面内の法線方向変化を二次項で補正し、元の頂点位置・法線を保ちつつ C¹ 連続な接触面を構成する手法である。
定式化
中間節点の生成： 両端頂点 x_i,x_j と中間節点を通る二次曲線の接線が、頂点法線 n_i, n_j と直交するように辺 i–j の中間節点 m_ij を決める。
\boldsymbol{m}_{ij} = \mathbf{C}_i^{(ij)}\,\boldsymbol{x}_i + \mathbf{C}_j^{(ij)}\,\boldsymbol{x}_j
本実装では、中間節点を法線方向にのみ移動する仮定のもと、正則化付き最小二乗法で補正係数行列を導出した。このとき辺の平均法線方向は m = (n_i + n_j)/‖n_i + n_j‖、a₁ = m·n_i、a₂ = m·n_j、曲率補正係数 K = 4/(16(a₁²+a₂²)+ε)（ε = 10⁻⁴）として：
\mathbf{C}_i^{(ij)} = \frac{1}{2}\mathbf{I} - K\,\boldsymbol{m}\otimes(a_2\boldsymbol{n}_j - a_1\boldsymbol{n}_i)
\mathbf{C}_j^{(ij)} = \frac{1}{2}\mathbf{I} + K\,\boldsymbol{m}\otimes(a_2\boldsymbol{n}_j - a_1\boldsymbol{n}_i)
補間行列 P： 中間節点は頂点座標の一次結合で表せるので、二次形状関数と組み合わせて、頂点-補間点の変換行列 P(ξ,η) を構成する。
\boldsymbol{r}(\xi,\eta) = P(\xi,\eta)\,\begin{bmatrix}\boldsymbol{x}_1\\ \vdots\\ \boldsymbol{x}_n\end{bmatrix}
入出力仕様
!CONTACT, INTERACTION=FSLID, GRPID=1, SMOOTHING=NAGATA, CONTACTPARAM=CPARAM1SMOOTHING に NAGATA を指定すると有効化される（デフォルト: NONE）。
実装の要点
- 要素計算側：中間節点を頂点座標の従属量として縮約した補間行列を使用（独立自由度の増加なし）
- 探索側：中間節点を明示的に生成し tri6/quad8 の形状関数で投影
- 法線は探索前に更新し MPI 集約で一貫性を確保
- Newton 反復中は幾何情報を固定（幾何揺れによる不安定化を防止）
精度向上結果
円筒-直方体転がり接触例題（SLAGRANGE, FSLID）で検証した結果をしめす。
メッシュ図：
転がり中の接触反力履歴。前半押し込み、後半が回転。SMOOTHINGありのほうが大幅に早く結果が収束する：
定量評価：
| ケース | 平均 Fy | 標準偏差 | 範囲 (max−min) | 基準解誤差 (%) |
|---|---|---|---|---|
| N=8 なし | 5.928 | 5.540 | 14.147 | 67.2% |
| N=8 NAGATA | 18.846 | 0.476 | 1.652 | 4.2% |
| N=16 なし | 11.470 | 1.360 | 4.051 | 36.6% |
| N=16 NAGATA | 18.158 | 0.110 | 0.383 | 0.4% |
| N=32 なし | 16.365 | 0.703 | 1.852 | 9.5% |
| N=32 NAGATA | 18.097 | 0.087 | 0.363 | 0.1% |
| N=64 なし | 17.617 | 0.095 | 0.461 | 2.6% |
| N=64 NAGATA | 18.069 | 0.096 | 0.355 | 0.1% |
| N=128 NAGATA（基準） | 18.081 | 0.086 | 0.358 | — |
N=16 で基準解誤差が 36.6% → 0.4%、標準偏差が 1.360 → 0.110 に改善。NAGATA なしでは N=64 でも 2.6% にとどまるのに対し、NAGATA 有では N=16 で 0.4% を達成しており、同等精度に必要な要素数を大幅に削減できる。
2. 接触反力ダンピング機能
問題設定
FREE/CONTACT の二値的な状態遷移は接触力の不連続な変化を招き、Newton 反復の収束性悪化や接触反力の振動につながる。
定式化
CONTACTNEAR 状態（投影点が存在し未貫通）にあるスレーブノードに対して、ダンピング残差・剛性を付加する。α = DAMP_ALPHA × r_ref、B_n は法線方向変位マッピングベクトル、Δu は増分変位ベクトルとして：
\Delta g_n = \boldsymbol{B}_n^T\,\Delta \boldsymbol{u}
\boldsymbol{r}_\mathrm{damp} = \alpha\, w(g)\, \Delta g_n\, \boldsymbol{B}_n
\mathbf{K}_\mathrm{damp} = \alpha\, w(g)\, \boldsymbol{B}_n \otimes \boldsymbol{B}_n
スムーズステップ重み関数 w(g)：
w(g) = \begin{cases} 1 & (g \le 0) \\ 1 - 3t^2 + 2t^3 & (0 < g < g_\mathrm{act}),\quad t = g/g_\mathrm{act} \\ 0 & (g \ge g_\mathrm{act}) \end{cases}
w(g) の C¹ 連続性により、接触状態変化時の残差の急激な変化を抑制する。
入出力仕様
!CONTACT, INTERACTION=FSLID, GRPID=1, SMOOTHING=NAGATA, DAMP_ALPHA=0.1, DAMP_GACT=0.01, CONTACTPARAM=CPARAM1| パラメータ | 意味 |
|---|---|
DAMP_ALPHA |
ダンピング係数（参照剛性に対する無次元比） |
DAMP_GACT |
ダンピング有効距離 |
実装の要点
- FREE/CONTACT → FREE/NEAR/CONTACT に状態拡張
- implicit/explicit 接触探索を統合（
scan_contact_stateにflag_ctAlgoをオプション引数で統合） - NEAR 状態では拘束を課さず、投影情報のみ保持してダンピング項を付加
- ダンピング項は既存の法線分配ベクトルを使用（幾何定義の二重化を回避）
動作確認結果
円筒-直方体転がり接触例題（N=16, NAGATA 有）で検証した結果：
| ケース | 平均 Fy | 標準偏差 | 範囲 (max−min) | 失敗数 |
|---|---|---|---|---|
| nodamp | 17.739 | 0.525 | 1.353 | 0 |
| damp_a01_g01 | 17.688 | 0.479 | 1.213 | 0 |
| damp_a10_g01 | 17.694 | 0.473 | 1.202 | 0 |
全ケースで発散なく収束し、振動範囲を最大 11%（1.353 → 1.202）低減。
3. MPI 並列動作確認
N=16, NAGATA 有, DAMP_ALPHA=0.1, DAMP_GACT=0.01 の条件で 1/2/4/8 MPI 並列を実行：
| ケース | 最大差 abs(Fy) | 平均差 abs(Fy) | ステップ数 |
|---|---|---|---|
| np1（基準） | — | — | 52 |
| np2 | 6.6e-6 | 1.3e-6 | 52 |
| np4 | 2.3e-6 | 4.2e-7 | 52 |
| np8 | 2.1e-6 | 4.2e-7 | 52 |
全並列数で同一ステップ数、機械精度レベル（最大差 6.6e-6）の数値一致を確認。
テスト
以下の単体テストケースを tests 配下に追加：
alag_nagata— Augmented Lagrange 法 + Nagata 静解析slag_nagata— Standard Lagrange 法 + Nagata 静解析slag_nagata_dynamic— Standard Lagrange 法 + Nagata 動解析
Edited by inagaki