ブシネスク問題（半無限弾性体の点荷重） — 理論と支配方程式

3Dソリッド要素の検証に非常に有用な古典的問題だ。弾性半無限体の表面に集中荷重 $P$ を加えたときの応力と変位の厳密解を1885年にBoussinesqが導出した。3D場の$1/r$特異性を含むから、要素の精度や特異点近傍のメッシュ設計を評価するのに最適だ。

主に2つ。ひとつはFEAの3Dソリッド要素のCode Verification。もうひとつはHertz接触解析の前段検証だ。Hertz接触理論はBoussinesq解の重ね合わせ（面圧分布の畳み込み積分）として導出されるから、Boussinesq問題をクリアできないソルバーで接触解析をやるのは危険だ。

荷重点を原点、荷重方向を$z$軸正方向とした円筒座標$(r, z)$で表す。$R = \sqrt{r^2 + z^2}$ として、変位場は

ここで $G = E/[2(1+\nu)]$ はせん断弾性係数。応力場で最も重要な$z$軸上の鉛直応力は

核心的な問いだ。FEAの離散化では特異点を正確に再現できない。だからこそ、荷重点から十分離れた位置で理論値との比較を行う。評価点は $z > 5h_{elem}$（$h_{elem}$は荷重点近傍の要素サイズ）を目安にする。荷重点での応力値はメッシュ依存で意味がないから評価対象外だ。

十分大きなモデル寸法を取ればよい。着目領域の最大距離を $L_{eval}$ とすると、モデルの外形は $R_{model} \geq 10 L_{eval}$、深さ $D_{model} \geq 10 L_{eval}$ にする。遠方の境界条件は底面を$z$方向固定、側面を半径方向ローラーにするのが標準だ。

軸対称性を活かしてCAX8要素（Abaqus）や2D軸対称要素で解くのが計算効率的だ。3D検証が必要なら1/4対称の3Dモデルを使う。

AbaqusのCINAX4（軸対称無限要素）やCIN3D8（3D無限要素）を外縁に配置すれば、モデルサイズを$R_{model} \geq 3 L_{eval}$程度まで縮小できる。計算コストが大幅に下がるから、パラメトリックスタディでは無限要素の活用を推奨する。

標準パラメータ: $P = 1$ MN、$E = 200$ GPa、$\nu = 0.3$。評価点は荷重点直下 $z = 0.1$ m。

メッシュサイズ5mmで1%以内の精度が得られる。3水準のメッシュでGCIを算出すると0.5%程度になり、十分に収束していることが確認できる。