ブシネスク問題（半無限弾性体の点荷重）

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}$ として、変位場は

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

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

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

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

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

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

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

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

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

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

以下の手順で進める。

評価点($z = 0.1$ m)での$\sigma_{zz}$と$u_z$が理論値の1%以内。メッシュ「細」での値を基準にGCI < 5%。この基準は片持ち梁ベンチマークと同じレベルだから、一貫性がある。

Hertz接触理論は、Boussinesq解の面圧分布に対する畳み込み積分として導出される。球と平面の接触で、接触面の半楕円面圧分布

その通り。V&Vの基本原則は「単純から複雑へ」だ。Boussinesq（線形弾性、軸対称）→Hertz（接触だが解析解あり）→一般接触（非線形、摩擦あり）と進む。各段階で参照解との一致を確認し、次の段階の信頼性基盤とする。

基礎設計の第一近似として使われる。地中の増加応力を計算し、圧密沈下量を推定する。具体的にはNewmarkチャートやFadum積分表を使って矩形分布荷重下の応力増分を求める手法が広く普及している。

作用応力が地盤の降伏応力に比べて十分小さい領域（常時荷重下の地盤）では弾性近似で実用的な精度が得られる。大変形や液状化を扱う場合は弾塑性モデルや有効応力解析に進む必要がある。その場合でもBoussinesq解は初期応力状態の推定やオーダー推定に使える。

軸対称CAX8Rモデルの場合の要点はこうだ。

OUTPUT, HISTORY, VARIABLE=PRESELECT で特定節点の変位を時刻歴出力。応力はEL PRINT で要素積分点値を取得するか、Python ODB APIで座標指定のフィールド出力を抽出する。probeValues メソッドで任意点の補間値を取得できるから、理論値との自動比較に便利だ。

CTRIAX6は6節点三角形で、四角形にする場合はTRAX3やTRAX6を使うが、Abaqusほど使いやすくはない。実務的には3D 1/4対称モデル（CHEXA-20）で解くことが多い。

Nastranの応力テンソルの出力順序は $\sigma_{xx}, \sigma_{yy}, \sigma_{zz}, \sigma_{xy}, \sigma_{yz}, \sigma_{xz}$ で、Abaqusは $\sigma_{11}, \sigma_{22}, \sigma_{33}, \sigma_{12}, \sigma_{13}, \sigma_{23}$ だ。特にせん断成分の順序が異なる（Nastranは $xy, yz, xz$、Abaqusは $12, 13, 23$）から、比較時に取り違えないように注意すること。

CalculiXはAbaqus互換の.inpファイルがほぼそのまま使える。C3D20要素で3D 1/4モデルを解く。注意点として、CalculiXのC3D20の積分点配置がAbaqusと微妙に異なることがある。結果は有効数字3〜4桁で一致するが、それ以上の精度が必要なら積分スキームの詳細をソースコード（src/e_c3d.f）で確認すべきだ。

ソースコードが公開されているため、要素定式化のアルゴリズムを第三者が直接検証できることだ。実際、CalculiXのC3D20の剛性マトリクス計算はe_c3d.f内のGauss積分で明示されている。商用ソルバーではこの透明性がない。学術論文でOSSベースの検証を行い、その結果を商用ソルバーの独立検証として引用するのは確立された手法だ。

Burmister（1945）が二層弾性体の解を導出している。Hankel変換を含む積分形式で、数値積分が必要だが閉じた形の参照解として使える。舗装設計（上層がアスファルト、下層が路盤）のFEA検証に使う。

完全接着条件なら共有節点で十分だが、すべり界面を考慮する場合は接触条件が必要になる。完全接着の検証ではBurmister解が参照になり、すべり界面の検証では別の解析解（層間摩擦なしのBurmister解）が必要だ。段階的に検証することが重要だ。

Lambの問題（1904）がそれに相当する。半無限体表面にステップ荷重を加えたときの過渡的弾性波伝播を扱い、P波・S波・Rayleigh波の到達時刻と振幅の理論解がある。

最短波長あたりのノード数が鍵だ。二次要素で1波長あたり5〜6要素、線形要素で10〜12要素が目安。S波の波長 $\lambda_s = V_s / f$ から必要な要素サイズが決まる。数値的な分散（周波数依存の位相速度ずれ）が許容範囲内かを確認するため、異なるメッシュ密度で波形を比較する必要がある。

弾塑性では厳密解が存在しないから、Boussinesqの弾性解をヤング率に関するパラメトリックスタディの基準として使う。荷重レベルが低い段階では弾性解に一致し、荷重増大に伴って弾塑性解が弾性解から乖離する。その乖離パターンが物理的に合理的か（圧縮降伏→塑性域拡大→支持力限界）を確認する。

まさにそうだ。単純な問題で基盤の信頼性を確立し、段階的に複雑化していく。Boussinesq → Hertz接触 → 多層地盤 → 弾塑性 → 動的問題という系統的なV&Vピラミッドの底辺にBoussinesqがある。この系統的アプローチこそがASME V&V 10とNAFEMS QSSが推奨する方法論だ。

固定境界からの反射で応力が過大評価される。底面固定が浅いと、荷重点直下の$\sigma_{zz}$が理論値より20%以上大きくなることもある。症状として、モデルサイズを変えると結果が変動する場合はこの影響を疑うべきだ。

無限要素の減衰関数が問題の特性と合っていない可能性がある。Abaqusの無限要素は指数減衰を仮定しているが、Boussinesq解は$1/r$減衰だから、厳密にはフィットしない。それでも実用上は十分な精度が得られるが、残差が気になる場合は有限要素領域を大きめに取る（無限要素との境界を遠ざける）のが安全だ。

それは正常な挙動だ。荷重点の応力はメッシュを細かくするほど際限なく増大する。これは物理的に正しい特異性の反映であり、FEAのバグではない。

低次要素（TET4、HEX8）で応力をコンタープロットすると、要素間で応力が不連続になりギザギザに見えることがある。これはC0連続の要素では応力（歪みの微分量）が要素境界で不連続になるためだ。

Zienkiewicz and Zhu（1992）が提案した事後誤差推定法だ。各節点のパッチ内の積分点応力から最小二乗法で節点応力を高精度に回復する。回復前後の差が局所的な離散化誤差の推定値になる。多くの商用ソルバーでオプションとして実装されている。Abaqusでは*ERROR INDICATOR で利用可能だ。