Claude CodeでCalculiXの入力ファイルを作り、片持ち梁を手計算で検算する

CalculiXの入力ファイル(.inp)は、Claude Codeに書かせること自体は難しくありません。難しいのは、解析が走って数字が出たあとに、その数字を信じてよいかを決めることです。CalculiXは単位を計算に使わないので、ヤング率の桁を間違えても、ソルバーは何も言わずに答えを返します。

この記事では、片持ち梁を題材に「手計算の目標値を先に固定してから、Claude Codeに入力ファイルを書かせ、出力と突き合わせる」流れを作ります。モデルは5要素のB32R梁で、手計算との差がどこまでなら正常かもマニュアルの数字で示します。

CalculiXの入力ファイルは何でできているか

CalculiXは有限要素法の構造解析プログラムで、GNU GPLで配布されています。ソルバーはAbaqus形式の入力を読むため、入力ファイルはAbaqusと同じ「* で始まるキーワードカード」の並びです。

マニュアルの構造は2段です。最初に節点・要素・材料・境界条件といったモデル定義、その後に1つ以上の *STEP ... *END STEP を置きます。モデル定義のカードは最初の *STEP より前に書きます。

*NODE PRINT と *EL PRINT は拡張子 .dat のASCIIファイルを作り、*NODE FILE と *EL FILE は可視化用の .frd を作ります。検算に使うのは .dat です。選択しなければ .dat は生成されません。

単位についてマニュアルの原則は1つです。入力する数字の単位が一貫していること。単位は計算に現れません。mm・N・MPaの組なら、長さはmm、力はN、応力とヤング率はMPa(N/mm²)で通します。

手計算の目標値を先に固定する

題材は、長さ100 mm、断面10 × 10 mmの鋼製片持ち梁です。固定端をz = 0に置き、自由端(z = 100)にx方向の集中荷重100 Nをかけます。ヤング率は210,000 MPa、ポアソン比は0.3です。

断面が正方形なので、曲げの向きによる断面二次モーメントの違いは考えなくて済みます。

• 断面二次モーメント: I = 10 × 10³ / 12 = 833.33 mm⁴
• 自由端のたわみ: δ = P L³ / (3 E I) = 100 × 100³ / (3 × 210000 × 833.33) = 0.1905 mm
• 固定端の最大曲げ応力: σ = M c / I = (100 × 100) × 5 / 833.33 = 60 MPa
• 固定端の反力と曲げモーメント: 100 Nと10,000 N·mm
• せん断応力の平均: P / A = 100 / 100 = 1 MPa

オイラー・ベルヌーイ梁の式にはせん断変形が入っていません。せん断補正係数を5/6、G = E / (2(1+ν)) = 80,769 MPaとして加えると、せん断分のたわみはPL / (κGA) ≈ 0.0015 mm、全体の約0.8%です。細長い梁では無視してよい大きさだと、この段階で見積もっておきます。

この値は、解析を回す前に check_beam.py に書き出してもらいます。

# check_beam.py - 理論値の単一の置き場(単位: mm, N, MPa)
L, b, h, P, E = 100.0, 10.0, 10.0, 100.0, 210000.0
I = b * h**3 / 12
 
tip_deflection = P * L**3 / (3 * E * I)   # 0.1905 mm
sigma_root_max = P * L * (h / 2) / I      # 60 MPa
 
def sigma_at(z, x):
    """自由端 z=L から測った曲げモーメントで求めた、位置 (z, x) の応力"""
    return P * (L - z) * x / I

数字を手で書かず関数にしておくのは、あとで積分点の位置に合わせた値を取り出すためです(次の節で使います)。

5要素のB32R梁をClaude Codeに書かせる

CalculiXの梁要素は3節点の2次要素で、内部で体積要素に展開されます。B32は20節点ブリック要素に、低減積分版のB32Rは低減積分ブリック要素に展開されます。マニュアルは、断面力が必要なときはB32を避けてB32Rを使うよう勧めています。

次のようなプロンプトで、手計算を渡したうえで入力ファイルを書かせます。

check_beam.py の値を目標に、CalculiXの入力ファイル beam.inp を書いて。
- 単位は mm, N, MPa。鋼 E=210000, nu=0.3
- z 方向 100 mm の梁を B32R 要素 5 本で。z=0 を全自由度固定
- z=100 の節点に x 方向 100 N
- 全節点の変位を *NODE PRINT で .dat に出す
- 固定端の節点集合 FIX の反力の合計も *NODE PRINT(TOTALS=ONLY, RF)で出す
書いたら ccx beam を実行し、自由端の ux を check_beam.py の
tip_deflection と比べて、差を % で報告して。

Claude Codeが書くデッキは、マニュアルの片持ち梁(B32R、*BEAM SECTION)と同じ骨格になります。5要素・11節点の場合は次の形です。

*NODE,NSET=NALL
1,0,0,0
2,0,0,10
3,0,0,20
4,0,0,30
5,0,0,40
6,0,0,50
7,0,0,60
8,0,0,70
9,0,0,80
10,0,0,90
11,0,0,100
*ELEMENT,TYPE=B32R,ELSET=EALL
1,1,2,3
2,3,4,5
3,5,6,7
4,7,8,9
5,9,10,11
*NSET,NSET=FIX
1
*BOUNDARY
1,1,6
*MATERIAL,NAME=STEEL
*ELASTIC
210000.,0.3
*BEAM SECTION,ELSET=EALL,MATERIAL=STEEL,SECTION=RECT
10.,10.
1.d0,0.d0,0.d0
*STEP
*STATIC
*CLOAD
11,1,100.
*NODE PRINT,NSET=NALL
U
*NODE PRINT,NSET=FIX,TOTALS=ONLY
RF
*EL PRINT,ELSET=EALL
S
*END STEP

*ELEMENT の節点並びは「端点・中間点・端点」です。要素2は節点3・4・5で、要素1の終点3を共有します。*BOUNDARY の 1,1,6 は節点1の自由度1〜6(並進3・回転3)の固定です。*BEAM SECTION の2行目は断面寸法、3行目は断面の向きを決めるベクトルで、マニュアルの例と同じ形に揃えています。

実行はジョブ名(拡張子なし)を渡します。

ccx beam
python3 check_beam.py

.dat の変位行は、節点番号に続けてvx・vy・vzが並ぶ形式です。マニュアルにある出力例は、displacements (vx,vy,vz) for set NALL and time 1. という見出しのあとに「節点番号vx vy vz」の行が続きます。この形式を前提に、比較のスクリプトをClaude Codeに足してもらいます。

# compare_dat.py - beam.dat の自由端変位を理論値と比べる
import re
from check_beam import tip_deflection
 
rows = {}
in_disp = False
for line in open("beam.dat"):
    if "displacements" in line:
        in_disp = True
        continue
    parts = line.split()
    if in_disp and len(parts) == 4:
        rows[int(parts[0])] = [float(v) for v in parts[1:]]
 
ux = rows[11][0]
err = (ux - tip_deflection) / tip_deflection * 100
print(f"ux(tip) = {ux:.5f} mm, 理論 {tip_deflection:.5f} mm, 差 {err:+.1f}%")

出力の書式はccxのバージョンや要素の種類で変わることがあるので、手元の beam.dat を一度開いて、見出しの文字列と列数が合っているかを確かめてから使います。合わなければ、スクリプトの側を直します。

差の許容幅はマニュアルの数字で決める

B32R梁の先端たわみは、手計算と完全には一致しません。体積要素の2次補間は2次の変形しか表せませんが、先端荷重の片持ち梁のたわみは3次の変化をするためです。

マニュアルの片持ち梁(100 mm、断面2 × 2 mm、先端に単位荷重)では次の差が出ています。

要素を増やすと誤差は縮みます。5要素の本稿のモデルでも、先端たわみは理論より数%小さく出るのが自然で、そこを「許容2〜3%」と置けます。逆に、差が1桁大きいときはモデルか入力の誤りを疑います。

反力を先に見る

たわみの前に、固定端の反力を見ます。釣り合いは要素の補間に左右されないので、反力は離散化の差がほとんど出ません。上のデッキでは *NODE PRINT,NSET=FIX,TOTALS=ONLY と RF で、固定端の節点集合の反力の合計だけを出しています。TOTALS=ONLYは個々の節点の寄与を省き、合計だけを印字する指定です。

マニュアルの別の片持ち梁の例では、.dat に次の形で出ます。

 total force (fx,fy,fz) for set FIX and time              0.1000000E+01

見出しの次の行に、fx・fy・fzの3つの値が並びます。この梁ならfxが荷重100 Nと釣り合う大きさになるはずで、fyとfzはほぼ0です。ここが100 Nから大きく外れていれば、たわみを見る前に荷重の節点・向き・固定の自由度を直します。

compare_dat.py には、次の数行を足してもらいます。

# 見出し行の次の空行を飛ばして、3つの値の行を読む
lines = open("beam.dat").read().splitlines()
i = next(k for k, l in enumerate(lines) if "total force" in l)
fx = float(lines[i + 2].split()[0])
print(f"RF fx = {fx:.3f} N, 理論 |100| N")

符号は荷重と逆向きになるのが自然ですが、書式と合わせて手元の beam.dat で確かめてください。*DLOAD のような分布荷重や *CLOAD を固定端の節点に掛けると、RFにその分が混ざります。今回は荷重を節点11だけに掛けているので、反力の合計がそのまま固定端の反力になります。

要素数を増やして差が縮むか確かめる

たわみの差が許容幅に入っていても、離散化による差かどうかは要素を増やさないと分かりません。5要素から10要素にして、差が縮むかを見ます。Claude Codeには次のように頼めます。

beam.inp を複製して beam10.inp を作って。
- 節点を 21 個(z=0,5,...,100)、B32R 要素を 10 本に変える
- 荷重は節点 21 の x 方向 100 N、固定は節点 1
- ccx beam10 を実行し、compare_dat.py と同じ比較で差を % で報告して
- 5要素と10要素の差を並べて、縮んでいるかを一行で書いて

マニュアルの表では、1要素の約10%が5要素で約2.4%に縮んでいます。10要素の差が5要素より小さければ離散化の差と読めます。変わらない、あるいは逆に広がるときは、メッシュではなく入力の問題です。このとき荷重の節点が自由端になっていない、といった節点番号の数え間違いが典型的な原因です。

応力は比較する場所をそろえる

応力の突き合わせでは、もう1つ落とし穴があります。*EL PRINT の S は積分点の応力で、固定端の最大応力60 MPaではありません。

低減積分の要素は、要素長の中でz方向に2点、断面内で ±1/√3の位置に積分点を持ちます。固定端の要素(z = 0〜20)なら、z = 4.226と15.774です。断面内のxは ±5/√3 = ±2.887 mm。

手計算の側も、その位置の値を出します。

from check_beam import sigma_at
print(sigma_at(z=4.226, x=5/3**0.5))   # 約 33.2 MPa

マニュアルの1要素・5要素の例は、この方法で積分点の応力を手計算と照合して、41.471 MPa(5要素)のように一致することを示しています。節点に外挿された値も .frd には入りますが、せん断応力は積分点では正しくても節点では外挿で外れた値になります。比べる位置を決めずに最大値同士を並べると、正常なモデルが不一致に見えます。

手計算が拾う桁の食い違い

手計算を先に持つ効果は、マニュアルの片持ち梁でも確かめられます。節の本文は「ヤング率100,000 MPa」と書き、表のたわみ2.50はP L³ / (3 E I) = 1 × 100³ / (3 × 100000 × 1.3333) = 2.50と一致します。一方、同じ節に掲載された入力デッキの *ELASTIC は 1E7,.3 です。

この値のまま解くと、たわみは100分の1の0.025付近になるはずです。本文と表の組が整合しているなら、デッキの桁が合っていません。ソルバーは単位を見ないので、この食い違いは手計算と並べて初めて現れます。

CLAUDE.mdと権限に書くこと

解析の検算を毎回頼めるようにするには、手順をプロジェクトのCLAUDE.mdに置きます。Claude Codeには、テストスイートやビルドの終了コード、出力を固定値と照らすスクリプトのような「読める検証手段」を渡すのが基本です。ここではテストの代わりに手計算スクリプトを置きます。

# CalculiX 解析のルール
- 単位は mm / N / MPa(= N/mm²)に固定。入力に単位は現れないので、
  材料・荷重・寸法はこの組で書く。
- .inp を書く前に、check_*.py に理論値を関数で書く。
- 解析後は ccx <job> → compare_dat.py を実行し、差を % で報告する。
- 差が 5% を超えたら、結果を直す前に入力(E の桁・荷重の向き・
  固定の自由度)を疑う。
- 応力は積分点位置で比べる。固定端の最大値と直接比較しない。

権限も、解析とスクリプトの実行に絞って許可します。

{
  "permissions": {
    "allow": [
      "Bash(ccx *)",
      "Bash(python3 check_beam.py)",
      "Bash(python3 compare_dat.py)"
    ]
  }
}

Bash(ccx *) は ccx で始まるコマンドに一致する形式の権限ルールです。権限ルールの書き方はClaude Codeのpermissionsコマンドの解説にあります。

同じ「移植先と元の数値を突き合わせる」型の検算は、MATLABからPythonへの移植でも使えます。形状側の寸法を確かめるなら、CadQueryの寸法アサーションが入口です。

結果が合わないときの切り分け

差が許容幅を超えたときは、次の順に確かめます。

実測の出力は、手元の ccx のバージョンと入力で変わります。この記事の数値は理論値とマニュアルの表に基づく目標で、実行結果そのものではありません。実行して出た数字を、上の許容幅に当てて判定してください。

まとめ

Claude Codeに解析を頼むときの要点は、入力ファイルの前に理論値を置くことです。ソルバーは単位を知らないので、答えの妥当性は手計算が担います。許容幅はマニュアルの数字(1要素で約10%、5要素で約2.4%)を目安にし、応力は積分点の位置で比べます。この3点をCLAUDE.mdに残しておけば、梁以外のモデルでも同じ検算の型を使い回せます。