ブシネスク問題（半無限弾性体の点荷重） — 数値解法と実装

$1/r^2$の応力勾配に対応するバイアスメッシュが必須だ。荷重点を中心とした放射状メッシュ（スパイダーウェブメッシュ）が最適で、幾何級数的に要素サイズを増大させる。

• 荷重点近傍: $h_{min} = L_{eval}/100$ 程度
• 外縁: $h_{max} = L_{eval}/2$ 程度
• バイアス比: 1.5〜2.0の幾何級数

軸対称の場合は荷重点（$r=0$軸上）に縮退した三角形要素を1列配置し、そこから四角形要素を放射状に展開する。

精度は出るが、同等精度に必要な要素数がHEX20の3〜5倍になる。検証目的ではマッピングメッシュ（六面体）を推奨する。ただし複雑形状への応用を見据えた自動テトラメッシュの精度検証をこの問題で行うのも有意義だ。

荷重点直下$z = 0.1$ mの$\sigma_{zz}$を指標に3水準で計算する。

観測収束次数: $p = \ln|(44.2-46.9)/(46.9-47.5)| / \ln 2 \approx 2.5$

GCI_fine: $1.25 \times |(-46.9-(-47.5))/(-47.5)| / (2^{2.5}-1) \approx 0.3\%$

特異点近傍では漸近領域に完全には入っていない段階で超収束が見られることがある。4水準目のメッシュ（$h_{min} = 2.5$ mm）を追加して$p$が安定しているか確認すべきだ。漸近領域での$p$安定が確認できれば、GCIの信頼性が保証される。

Abaqus: CAX8Rが定番。*CLOAD で軸上節点に荷重印加。底面の無限要素はCINAX4。後処理はPythonスクリプトでODB APIから応力抽出可能。

Nastran: CTRIAX6で軸対称解析可能だが使いにくい。3DモデルでCHEXA-20が実用的。SOL 101、GPSTRESS出力で任意点の応力テンソルを取得。

CalculiX: Abaqus互換のC3D20を使用。軸対称要素はないから3Dモデルで解く。ccx で実行し、.frdをParaViewで可視化。

Code_Aster: MODELISATION='AXIS'で軸対称。QUAD8相当の要素。SALOMEでメッシュ生成、コマンドファイルで境界条件設定。

同等のメッシュ密度・要素次数なら変位で4〜5桁、応力で3〜4桁一致する。差異は主に節点外挿アルゴリズムの違いに起因する。積分点での値で比較すれば一致度はさらに向上する。