無料のCAEは昔からあります。CalculiXもOpenFOAMも、実務に十分な力を持っています。それでも現場で使われにくい理由は、ソフトの性能よりも、出てきた数字を信じていいのか、自分で判断しにくいことだと考えています。
NovaSolverでは、オープンソースのソルバを束ねた解析ソフト NovaSolver Studio を作っています。この記事では、精度の話だけをします。解析解とどれだけ合うかと、合わないときに何を言ってくるかです。
何を解けるか
| 解析 | ソルバ | 確かめ方 |
|---|---|---|
| 構造（線形静解析・固有振動・熱応力・自重/遠心力） | 内蔵FEM（2次四面体）／CalculiX | はり理論、Heywood の Kt、メッシュ収束（GCI） |
| 構造（接触・大変形・弾塑性・ゴム） | 内蔵FEM | Lamé、Euler エラスティカ、塑性崩壊荷重 |
| 熱（定常伝導） | 内蔵FEM／CalculiX | フーリエの法則、フィン理論、GCI |
| 電気（静電場・直流・ジュール発熱） | 内蔵FEM | 平行平板・同軸の容量、両端冷却棒の温度上昇 |
| 内部流れ | OpenFOAM | Ghia ら（キャビティ）、Poiseuille、質量収支 |
| 外部流れ | OpenFOAM（snappyHexMesh＋壁面層） | 球の抗力（Schiller–Naumann と +0.3%） |
| 水冷（流体と固体の連成伝熱） | OpenFOAM（chtMultiRegionSimpleFoam） | 熱収支 |
形状は STEP／IGES／BREP／STL を読み込めます。
結果には、何を確かめたかがそのまま書かれた判定の行がつきます。
ran=yes / converged=yes / mesh convergence not checked
- 解析解がある例題では、その解析解と照合し、合えば
matches referenceと返します。 - 任意形状の部品では、粗い・標準・細かいの3つのメッシュで解き直し、追跡する量（構造では最大変位、外部流れでは抗力係数）について Roache の GCI を求めます。GCI ≤ 5%、全段でソルバが収束、段ごとに要素数が1.3倍以上、値の変化が落ち着く傾向にある、をすべて満たすと
mesh convergence (GCI) passed、満たさなければmesh convergence not reachedです。
以前はこれを「検証済み」と1語で返していましたが、やめました。メッシュ収束を確かめても分かるのは対象量のメッシュによる誤差の見積もりで、GCI も真の誤差の上限ではありません。材料・荷重・拘束のモデル化が妥当か、実物と合うか、安全か、までは分からないためです。
21の解析解・文献値との比較
冒頭の図が、配布物に同梱している21の検証例の結果です。参照値は解析解のほか、Ghia らの文献値や Heywood の近似式を含みます（Studio 0.6.2）。灰色の帯は、問題ごとに決めた許容範囲です。
全問が許容範囲に入りました。最大は、はりの塑性崩壊荷重の +2.44% です。参照値は理想剛塑性の崩壊荷重 P = σy·b·h²/(4L) で、Studio の値は最後に収束した荷重増分です。正の誤差なので耐荷力を少し過大に見積もっています。メッシュと荷重増分のどちらの影響が大きいかは、まだ切り分けていません。なお、この例は崩壊点まで載荷するため、判定上は「収束」せずに終わる設計です。
誤差の大きい順に並べると、次のとおりです。
| 問題 | 参照値 | Studio | 差 |
|---|---|---|---|
| はりの塑性崩壊荷重 | 1250 N | 1280.5 N | +2.44% |
| 片持ち梁の先端たわみ（せん断補正つき） | 0.9831 mm | 0.9755 mm | −0.77% |
| 穴あき板の Kt（Heywood） | 2.512 | 2.531 | +0.75% |
| 大変形片持ちはり（PL²/EI = 1） | 60.34 mm | 60.22 mm | −0.20% |
| 回転円板の中心応力 | 35.51 MPa | 35.58 MPa | +0.20% |
| 圧入（Lamé） | 76.88 MPa | 76.96 MPa | +0.11% |
キャビティ流れ（Ghia ら, Re=100）は、中心線上の u 分布の平均絶対誤差がふた速度基準で 0.12% でした（中心線上の最小 u は −0.2109 に対して −0.2115 m/s）。
開発の途中では、NAFEMS のベンチマークと ISO 10211 の熱橋も解いています（0.2.2開発版）。
- NAFEMS LE1：参照解との差 0.4% 以内
- NAFEMS LE10：0.2% 以内
- NAFEMS T4：温度差 約0.05 K
- ISO 10211 の2問：規定（温度差0.1 K以内、熱流差の範囲）を満たした
合わないときに何を言うか
解析解がある問題で合うのは、ある意味で当然です。実務で困るのは解析解がない部品で、そこでは「収束しない数字」をどう扱うかが問題になります。
角を丸めていないL字の金具（アルミ、壁側の面を固定、穴に下向き500 N）を、3段のメッシュで解かせると、Studio 0.7.5 は次のように返しました。
ran=yes / converged=yes / mesh convergence not reached
最大ミーゼス応力: 61.77 MPa
最大変位: 0.3145 mm
WARNING: 最大応力はメッシュを細かくするほど上がり続けています（49 → 50 → 62 MPa）。
鋭い角や完全固定した面の縁など、理想化した形・拘束による応力特異性の可能性があり、
この最大値そのものは強度の判断に使えません。
この警告は、3段のメッシュでピーク応力が増え続け、最後の増分が3%を超えたときに出る検出則です。特異性を数学的に証明するものではなく、「疑わしい」という合図です。
角と穴のまわりだけメッシュを細かくして解き直すと、最大変位の GCI は 0.33% になり、mesh convergence (GCI) passed になりました。一方で角の最大応力は 63 → 77 → 94 MPa と増え続け、警告は消えません。
| 量 | 値 | 判断 |
|---|---|---|
| 最大変位 | 0.318 mm | GCI 0.33%（この量のメッシュ収束は基準を満たした） |
| 角の最大応力 | 63 → 77 → 94 MPa | 特異性が疑われる。この値では判断しない |
| 穴の面の最大応力 | 16.645 → 16.653 MPa | 局所細分化の前後で 0.05% の変化（3段の収束評価は未実施） |
鋭い再入隅角や拘束端の応力集中は、教科書に書いてあっても、締め切りに追われていると見落とします。その兆候をソフトの側で知らせるのが狙いです。
実務寄りの3題
水冷コールドプレート：熱収支で確かめる
GPU 700 W を、銅板の中の6本の水路（3 × 4 mm）で冷やす形です。台座面に700 Wを一様にかけ、水と銅の連成伝熱を約26万要素で解きました。GPU内部やTIM、接触熱抵抗は含めていません。
| 水量 | 銅板の最高温度 | 出口水温 | 水が運んだ熱 | 手計算の出口水温 |
|---|---|---|---|---|
| 0.8 L/min | 90.7 ℃ | 37.3 ℃ | 681.8 W（97.4%） | 37.6 ℃ |
| 0.4 L/min | 103.5 ℃ | 49.9 ℃ | 692.8 W（99.0%） | 50.2 ℃ |
出口水温は Q = ṁ c ΔT の手計算と約0.3 Kの差、熱収支の未閉合は2.6%／1.0%です。ただし、熱収支が閉じていることは整合性の確認にすぎません。壁面熱伝達や銅板の最高温度の精度を、独立に保証するものではありません。
水量を半分にすると、水の温度上昇はほぼ2倍（12 → 25 K）になります。一方、銅板の最高温度は12.8 Kしか上がりません。また、水路ごとの流速には明らかな偏りがあり、入口側の水路は出口側よりかなり遅くなっています。分配溝の形を変えた比較はまだ解いていないので、改善の候補として次に確かめるところです。
EVバスバー：ジュール発熱は I² に比例するか
銅バスバーとアルミ端子を M6 ボルト1本で留めました。接触面に0.2 mmの等価層を置き、締結部の接触抵抗として良好 10 µΩ／緩み 200 µΩ を与えています（典型値の仮定で、締付け力から求めた値ではありません）。物性と熱伝達率は温度によらず一定です。
全体の発熱は手計算の I²R と一致し、温度上昇の比は 300 A／150 A で 3.99 でした。この線形定常モデルでは理論どおり約4倍です。実物では抵抗の温度依存や接触状態の変化が加わります。
鋼製ラックの固有振動数
幅900 × 奥行450 × 高さ1800 mm の一体のラックで、脚の底を完全固定しました。1次の固有振動数は 9.73 Hz、さらに上面の水平変位を拘束すると 36.77 Hz でした（柱は中実の40 mm角）。振動数の比較だけで、地震応答や転倒は解いていません。
160年前の特許を解く
参照解のない問題でどこまで言えるかを試す題材として、パブリックドメインになっている古い米国特許を選びました。
大砲つき鋤（US 35,600、1862年）
寸法は図面の比率から仮定し、薬室に 60 MPa を静的にかけ、刃を完全固定しています。
- 砲身の周方向応力は約79 MPaでした。Lamé の厚肉円筒の式で見積もった約81 MPa（内半径25 mm、外半径65 mm、内面）と同程度です。
- 支持反力 118 kN は、ガス圧 × 口径面積の 117.8 kN と一致します。
- 画面の最大値 687 MPa は、刃を完全固定した辺に出るピークで、判定には使えません。
一瞬の衝撃を静的にかけた計算で、鋤の慣性や土の抵抗は入っていません。
エンジンのない飛行船エアリオン（US 43,449、1864年）
特許には「浮力の余り1ポンドで時速1マイル、200ポンドなら時速200マイル」と書かれています。機首10°・5 m/s の1条件で外部流れを k-ω SST で解き、Cd = 0.238、Cl = 0.563 を求めました。この係数を一定として、浮力の余りと空気力を釣り合わせ、水平速度を推定しています。
この簡略モデルでは、1ポンドで時速約0.9マイル、200ポンドで約12.5マイルです。空気力が速度の2乗に比例するなら、速度は力の平方根でしか増えません。200倍の力でも約14倍です。
ただし、係数を求めた5 m/sと、1ポンドのときの速度（約0.4 m/s）は大きく離れています。ゴンドラや網の抵抗も省いています。今回は1つのメッシュで解いただけで、参照値もメッシュ収束の確認もないため、Studio の判定は「実行・収束」までです（外部流れでも、球の抗力の相関式との照合や、抗力係数のメッシュ収束を確かめる経路はあります）。
使い方
- GitHub Releases から zip を落として、
NovaSolverStudio.exeを起動します。 - Python があれば
uvxでも動きます（README 参照）。 - MCP に対応しているので、Claude などから「この STEP を読み込んで底を固定、上から 500 N。メッシュ収束まで確かめて」と頼めます。
ライセンスは GPL-3.0 です。こういう形の検証例がほしい、という要望や、合わなかった例の報告を歓迎します。