摩擦判定における最大静止摩擦力を拡張ステップ内で固定する現行仕様だと、摩擦係数が荷重経路の勾配を超えるとALAGの拡張ループが収束しない (!615 の議論) ## 用語 - **拡張ループ**:拡張ラグランジュ法における，乗数λを更新して Newton を解き直す外側の反復。`AUGITER` 回まわす。 - **摩擦力の上限 (摩擦錐の半径)**: 最大静止摩擦力=摩擦係数 × 法線接触力。すべり状態節点の摩擦力はこの値になる。 - **現行仕様 (FROZEN)**: この上限を計算する法線接触力に「前回の拡張で確定した乗数」を使う。拡張ステップの途中では上限を動かさない。 - **新仕様 (FOLLOW)**: 法線接触力に「今この Newton 反復で接触面に実際に働いている力 (乗数 + ペナルティ × 貫入)」を使う。上限が貫入に応じて動く。 - **荷重経路の勾配**: 外力を接触面が受け止めるとき、「面を押し付ける成分」が「面に沿う成分」に対してどれだけあるかの比。斜面を水平力で押し上げる楔 (下図) では、水平力 F が斜面に対して 法線方向 F·sinθ と 接線方向 F·cosθ に分かれるので、比は tanθ = 斜面の勾配そのもの。この比が小さいほど荷重の大半を摩擦力 (面に沿う力) が支えることになり、摩擦力の上限の決め方が結果を左右する。本文の「摩擦係数 ÷ 荷重経路の勾配」は、摩擦が支える分と押し付けが支える分の比。 {width=900 height=505} --- ## 0. 要約 摩擦つき点-面の拡張ラグランジュで、荷重を接触が斜めに支える経路 (斜面を押し上げる楔など) では、現行仕様は拡張ループが 2 つの状態を往復して収束しない。これは実装の誤りではなく、拡張ステップ内で摩擦力の上限を固定する設計選択の帰結で、「摩擦係数 ÷ 荷重経路の勾配 > 1」で起きる。上限をその時点の押し付け力から決める新仕様 (!615) にすると同じ deck が解け、既存の in-tree deck も悪化しない。既定をどちらにするかを相談したい。 ## 1. 現状の問題 - **何が起きるか**: 楔 (斜面勾配 0.115、摩擦係数 0.2、595 節点) を力制御で押し上げると、各接触反復の `AUGITER` 回の拡張のあいだ乗数が 2 つの状態を往復し続けて確定せず、接触反復が収束しないまま停止する。\ tests 向けの縮小版 (dt 0.1 × 10) では t = 0.1 で "Contact failed to Converge"。摩擦係数 0.12 / 0.2 / 0.3、固定刻みと自動刻みの楔 7 種類はいずれも完走せず、停止時点の合計接触法線力は解析つり合い (5624.1 N) の −17〜−40%。 - **どんな条件で起きるか**: 摩擦係数 ÷ 荷重経路の勾配 が 1 を超えるとき。既存の in-tree deck は変位制御で、この比が小さいので表に出なかった (摩擦係数を 1.2 に上げても完走する)。 - **再現**: !615 に含めた縮小楔 `contact_alag_iter/wedge_conefollow_iter` (5 秒) を、`FRICTION_CONE=FOLLOW` を付けて／外して回すとその場で比べられる。 ## 2. 現行仕様の意図 - 拡張ステップ内で摩擦力の上限を固定すると、接線剛性が対称に保てる (Laursen–Simo の algorithmic symmetrization と同型)。Coulomb 則との整合は拡張ループに任せる設計 (e9e0208a、cde14769、!613)。 - 利点: 対称 solver で済み、実装が単純。 - 限界: 上限が押し付け力の変化に 1 拡張ぶん遅れる。斜面勾配 s の荷重経路では、1 節点のつり合いを書くと乗数の誤差が拡張 1 回ごとに −(摩擦係数 / s) 倍される (詳細は付録)。摩擦係数 < s なら符号を変えながら縮み、摩擦係数 > s なら発散して 2 状態の往復に落ちる。勾配 0.06〜0.30 × 摩擦係数 0.05〜0.4 の楔 9 本で、この倍率の予測と実測の収縮率が 2〜3 桁一致した (例: 予測 0.870 / 実測 0.869)。ペナルティ係数には依存しない。 ## 3. 新仕様と、その利点・根拠 - **何を変えるか** (3 点をセットで): 1. 摩擦力の上限を、その Newton 反復で実際に働いている押し付け力 max(0, 乗数 + ペナルティ × 貫入) から決める。 2. 上限が貫入で変わるぶんの線形化 (滑り摩擦力の結合ブロック) を接線剛性に足す。この項の「滑り方向」には、反復ごとに振れる試行摩擦力の向きでなく、拡張ステップ内で固定されている接線乗数の向きを使う\ (対称化と同じ「固定量で線形化する」思想を向きにだけ適用。収束点では両者が一致する)。 3. 結合ブロックは非対称なので、摩擦つき拡張ラグランジュの線形 solver を一般行列で初期化する。標準ラグランジュが摩擦つきで既に使っている判定と同じ経路。\ 半分だけ (上限の追従だけ、結合ブロックだけ) では Newton が壊れるので、1 つのキーワードでまとめて切り替える。 - **根拠** (数字は !615 の検証表と手元の台帳): | 検証モデル | 現行 (FROZEN) | 新 (FOLLOW) | |---|---|---| | 楔 (勾配 0.115、摩擦 0.2、40 substep) | 途中で停止、真値 −30% | 完走、解析つり合いと 0.001% | | 楔 縮小版 (tests、10 substep) | t = 0.1 で失敗 | 完走、0.018% | | 固着 patch test (拡張 20 回) | R_x 13.3145 N | 同じ値 (相対差 1.6e-8) | | in-tree 3beam_fslid、摩擦 1.2 | 完走、Newton 113 回 | 完走、Newton 82 回 | | in-tree alag_nagata (平滑化) | 完走、Newton 42 回 | 完走、Newton 52 回 (cutback なし) | | 実機規模のモデル (数万節点、圧入、摩擦 0.12) | 完走 (減衰などの対処つき) | 完走、Newton 約半分、解の差 0.03% 以内 | | 摩擦なしの deck 9 件 | — | 現行と bit 一致 | - **実装の形** (現在の !615): `!CONTACT_ALGO, TYPE=ALAGRANGE, FRICTION_CONE=FROZEN|FOLLOW`、既定は FROZEN で既存の結果は bit 不変。 ## 4. 相談したいこと 1. **既定を FOLLOW にしてよいか**。手元の in-tree deck では FOLLOW が劣る例は見つかっていない (○× 表を後述)。代償は非対称 solver で 1 反復あたり約 +27% の時間と、CG 系の deck では BiCGSTAB 等への変更が要ること。 2. **現行仕様 (FROZEN) を残す必要があるか**。上位互換なら残さない選択もある。残すなら opt-out の形。 3. 既定を変える場合の **参照解の扱い** (摩擦つき deck の参照解が動く)。 ○× 表 (in-tree の拡張ラグランジュ deck 全 16 件 + 摩擦係数 1.2 の変種 3 件。摩擦つきの CG 指定 deck は両仕様とも BiCGSTAB で比較。Newton は run 全体の合計回数。「収束解」= AUGITER=50 で拡張ループを回し切った解で、その行の「差」は両仕様の収束解どうしの最大変位の相対差): **摩擦つき deck (deck 既定の拡張回数 と AUGITER=50 の 2 行ずつ。右の 3 列は最終時刻の最大変位 [mm] を 現行 / 新 / 差 (現行 − 新) で書く。摩擦つきの CG 指定 deck は両仕様とも BiCGSTAB で比較)** | deck | 摩擦係数 | AUGITER | 現行 FROZEN: 完走 / Newton | 新 FOLLOW: 完走 / Newton | 最大変位 現行 | 最大変位 新 | 差 (現行 − 新) | |---|---|---|---|---|---|---|---| | 2box_hex_fric_iter | 0.3 | 3 (deck 既定) | ○ / 9 | ○ / 8 | 0.1018566 | 0.1018035 | +5.3e-5 | | | | 50 | ○ / 59 | ○ / 56 | 0.1018138 | 0.1018138 | +9e-15 | | 2box_hex_fric_iter | 1.2 | 3 | ○ / 9 | ○ / 8 | 0.1018187 | 0.1018035 | +1.5e-5 | | | | 50 | ○ / 58 | ○ / 56 | 0.1018138 | 0.1018138 | +1e-16 | | 3beam_fslid_iter | 1.2 | 3 | ○ / 113 | ○ / 82 | 0.2183696 | 0.2182989 | +7.1e-5 | | | | 50 | ○ / 1306 | ○ / 1241 | 0.2182171 | 0.2182169 | +1.8e-7 | | 3beam_fslid_iter | 0.2 | 3 | ○ / 84 | ○ / 82 | 0.2246201 | 0.2246209 | −8e-7 | | alag_nagata (平滑化) | 0.2 | 2 (既定) | ○ / 42 | ○ / 52 | 0.9912141 | 0.9912141 | +1.6e-8 (節点ごとの差の最大は 1.4e-3 相対) | | | | 50 | ○ / 803 | ○ / 960 | 0.9912140 | 0.9912141 | −1.1e-7 (節点ごとの差の最大は 7.9e-4 相対、50 回でも残る) | | 2beam_fslid_iter | 0.2 | 3 | ○ / 56 | ○ / 62 | 3.062773 | 3.062773 | +2.8e-7 | | | | 50 | **✗ t=0.72** / 7485 | **✗ t=0.78** / 10878 | — | — | (両仕様とも途中で止まる) | | 2beam_fslid_iter | 1.2 | 3 | ○ / 68 | ○ / 68 | 3.061421 | 3.061391 | +3.1e-5 | | | | 50 | **✗ t=0.72** / 7729 | **✗ t=0.85** / 11065 | — | — | (同上) | | 2beam_sslid_iter | 0.2 | 3 | ○ / 38 | ○ / 36 | 0.5000000 | 0.5000000 | +2e-9 | | cgs4_alag | 0.3 | 既定 | ○ / 294 | ○ / 261 | (節点ごとの差の最大 2.5e-5 相対) | | | | wedge_conefollow_iter (新規) | 0.2 | 既定 | **✗ t=0.1** / 230 | ○ / 176 | — | (解析つり合いとの差 1.8e-4) | — | **摩擦なし 5 件 + TIED 4 件**: 両仕様で結果が bit 一致 (FOLLOW は摩擦ゼロでは何もしない)。 読み方: - FOLLOW で完走しなくなった deck は無い。Newton 回数は 0.73〜1.24 倍。 - 既定の拡張回数で最大変位に差がある deck (2box、3beam@1.2) は、AUGITER=50 では差が 1e-7 以下に消える。 - alag_nagata は最大変位が強制変位で決まるので最大値の差は小さいが、節点ごとに見ると 50 回でも 8e-4 (相対) の差が残る (現行側が収束しない台)。2beam_fslid は 50 回だと両仕様とも途中で止まる。この 2 台は優劣を判定できない。 - **上位互換にならない条件が 1 つ**: solver に CG を指定した摩擦つき deck を deck のまま FOLLOW にすると、非対称行列で CG が停滞して 8 件中 4 件が異常終了・2 件が途中停止。既定を FOLLOW にするなら、摩擦つき拡張ラグランジュでは solver を BiCGSTAB 等へ切り替える必要がある (deck 側の変更か、コード側で自動選択・警告か)。