GitHub Copilot (Claude)でCFDのコード生成
低マッハ数近似というOpenFOAMの伝熱系のソルバとかの基本になるモデルがあります。それ使ったソルバをGitHub Copilot上で生成AI（Claude Sonnet）を動かして，作ってみました。こっちは指示しか与えませんでした。パラメータの設定とかは修正しましたが基本ちゃんと動きます。ブジネスク近似が破綻するような温度差でも対応可。
高温度差の熱対流の解析コード （ここをクリックするとソースファイルの詳細が出ます）
"""
2次元レイリーベナール対流解析コード（Julia版・次元あり）
- 低マッハ数近似（Low-Mach Number Approximation）
- PISO法（Pressure Implicit with Splitting of Operators）
- スタガード格子
- 下壁：高温（T=350K），上壁：低温（T=300K）
- 全壁：すべりなし境界条件
- 作動流体：空気（理想気体）
【無次元版からの主な変更点】
1. 全物理量を次元あり（SI単位系）で扱う
2. 密度：ρ = p0 / (R * T)（理想気体，p0 は熱力学的圧力）
3. 粘性係数 μ，熱拡散係数 α は定数（参照温度での値）
4. 浮力項：(ρ - ρ_ref) * g
5. 圧力ポアソン：∇·(1/ρ ∇p') = (∇·u* - S_T) / dt
6. 密度変化ソース項：S_T = α ∇²T / T
"""
using Printf
using LinearAlgebra
using Statistics
# ============================================================
# 物理定数・流体物性（空気，参照温度 T_ref = 325 K）
# ============================================================
const g = 9.81 # 重力加速度 [m/s²]
const R_gas = 287.0 # 空気の気体定数 [J/(kg·K)]
const c_p = 1005.0 # 定圧比熱 [J/(kg·K)]
# 温度条件
const T_hot = 450.0 # 下壁温度 [K]
const T_cold = 300.0 # 上壁温度 [K]
const T_init = 300.0 # 初期温度 [K]
const T_ref = 0.5 * (T_hot + T_cold) # 参照温度 [K] = 325 K
# 参照状態（T_ref における空気物性）
const p0 = 101325.0 # 熱力学的圧力（大気圧，空間一様）[Pa]
const rho_ref = p0 / (R_gas * T_ref) # 参照密度 [kg/m³] ≈ 1.084
const mu = 1.96e-5 # 動粘性率（粘性係数）[Pa·s] （325K の空気）
const nu = 1e-2 #mu / rho_ref # 動粘性係数 [m²/s]
const lambda = 0.0281 # 熱伝導率 [W/(m·K)] （325K の空気）
const alpha = 1.e-2 #lambda / (rho_ref * c_p) # 熱拡散係数 [m²/s]
# 粘性係数が大きすぎると自然対流が起きない
# ============================================================
# 計算領域・格子
# ============================================================
const Lx = 2.0 # x方向長さ [m]
const Ly = 1.0 # y方向長さ [m]（高さ）
const Nx = 64 # x方向セル数
const Ny = 32 # y方向セル数
const dx = Lx / Nx # [m]
const dy = Ly / Ny # [m]
const xc = [(i + 0.5) * dx for i in 0:Nx-1] # セル中心 x座標 [m]
const yc = [(j + 0.5) * dy for j in 0:Ny-1] # セル中心 y座標 [m]
# ============================================================
# 時間・ソルバパラメータ
# ============================================================
const dt = 5e-4 # 時間刻み [s]
const t_end = 10 #100.0 # 終了時間 [s]
const save_interval = 1000 # 結果保存間隔（ステップ数）
const n_correctors = 2 # PISO 圧力修正回数
const n_pressure_iter = 1500 # 圧力ポアソン反復数
const pressure_tol = 0.05 # 圧力収束判定 [Pa]
const sor_omega = 1.6 # SOR 緩和係数
# ============================================================
# 参照パラメータの表示用（レイリー数・プラントル数）
# ============================================================
const beta = 1.0 / T_ref # 体積膨張係数（理想気体）[1/K]
const Ra = g * beta * (T_hot - T_cold) * Ly^3 / (nu * alpha)
const Pr = nu / alpha
# ============================================================
# 密度の計算
# 低マッハ数近似：ρ = p0 / (R * T)
# ============================================================
"""セル中心密度を計算（低マッハ数近似：ρ = p0 / (R_gas * T)）[kg/m³]"""
function compute_density!(rho, T)
@inbounds for j in 1:Ny, i in 1:Nx
rho[i, j] = p0 / (R_gas * max(T[i, j], 1.0)) # ゼロ除算防止
end
end
"""u面上の密度を線形補間 [kg/m³]"""
function interpolate_rho_u(rho)
rho_u = fill(rho_ref, Nx+1, Ny)
@inbounds for j in 1:Ny
rho_u[1, j] = rho[1, j]
rho_u[Nx+1, j] = rho[Nx, j]
for i in 2:Nx
rho_u[i, j] = 0.5 * (rho[i-1, j] + rho[i, j])
end
end
return rho_u
end
"""v面上の密度を線形補間 [kg/m³]"""
function interpolate_rho_v(rho)
rho_v = fill(rho_ref, Nx, Ny+1)
@inbounds for i in 1:Nx
rho_v[i, 1] = rho[i, 1]
rho_v[i, Ny+1] = rho[i, Ny]
for j in 2:Ny
rho_v[i, j] = 0.5 * (rho[i, j-1] + rho[i, j])
end
end
return rho_v
end
# ============================================================
# 境界条件
# ============================================================
"""すべりなし・不透過境界条件"""
function apply_velocity_bc!(u, v)
u[1, :] .= 0.0
u[Nx+1, :] .= 0.0
v[:, 1] .= 0.0
v[:, Ny+1] .= 0.0
end
"""温度のゴーストセルを返す [K]
下壁：T = T_hot（ディリクレ）
上壁：T = T_cold（ディリクレ）
左右壁：断熱（ノイマン）
"""
function get_T_extended(T)
T_ext = zeros(Nx+2, Ny+2)
T_ext[2:Nx+1, 2:Ny+1] .= T
T_ext[2:Nx+1, 1] .= 2.0 * T_hot .- T[:, 1] # 下壁
T_ext[2:Nx+1, Ny+2] .= 2.0 * T_cold .- T[:, Ny] # 上壁
T_ext[1, 2:Ny+1] .= T[1, :] # 左壁（断熱）
T_ext[Nx+2, 2:Ny+1] .= T[Nx, :] # 右壁（断熱）
T_ext[1, 1] = T[1, 1]
T_ext[Nx+2, 1] = T[Nx, 1]
T_ext[1, Ny+2] = T[1, Ny]
T_ext[Nx+2, Ny+2] = T[Nx, Ny]
return T_ext
end
# ============================================================
# 中間速度の計算
#
# u成分：ρ Du/Dt = -∂p/∂x + μ ∇²u
# v成分：ρ Dv/Dt = -∂p/∂y + μ ∇²v + (ρ - ρ_ref) * g ← 浮力
#
# 離散化：u* = u + (dt / ρ) * RHS
# ============================================================
"""u速度の中間値（u*）を計算 [m/s]"""
function compute_u_star!(u_star, u, v, p, rho_u)
fill!(u_star, 0.0)
@inbounds for j in 1:Ny, i in 2:Nx
u_c = u[i, j]
rho_f = rho_u[i, j] # [kg/m³]
# 対流項 [m/s²]
u_e = (i < Nx) ? 0.5*(u[i, j] + u[i+1, j]) : 0.0
u_w = (i > 2) ? 0.5*(u[i-1, j] + u[i, j]) : 0.0
conv_x = (u_e^2 - u_w^2) / dx
if j < Ny
v_n = 0.5*(v[i-1, j+1] + v[i, j+1])
u_n = 0.5*(u[i, j] + u[i, j+1])
else
v_n = 0.0; u_n = 0.0
end
if j > 1
v_s = 0.5*(v[i-1, j] + v[i, j])
u_s = 0.5*(u[i, j-1] + u[i, j])
else
v_s = 0.0; u_s = 0.0
end
conv_y = (v_n*u_n - v_s*u_s) / dy
# 拡散項（μ/ρ = ν）[m/s²]
diff_x = (u[i+1, j] - 2*u[i, j] + u[i-1, j]) / dx^2
u_jm1 = (j > 1) ? u[i, j-1] : -u[i, j]
u_jp1 = (j < Ny) ? u[i, j+1] : -u[i, j]
diff_y = (u_jp1 - 2*u[i, j] + u_jm1) / dy^2
# 圧力勾配項（1/ρ * ∂p/∂x）[m/s²]
dp_dx = (p[i, j] - p[i-1, j]) / dx # [Pa/m]
u_star[i, j] = u_c + dt * (
-(conv_x + conv_y)
+ nu * (diff_x + diff_y)
- dp_dx / rho_f
)
end
u_star[1, :] .= 0.0
u_star[Nx+1, :] .= 0.0
end
"""v速度の中間値（v*）を計算 [m/s]"""
function compute_v_star!(v_star, u, v, p, T, rho_v)
fill!(v_star, 0.0)
# 温度をv面へ線形補間 [K]
T_v = zeros(Nx, Ny+1)
@inbounds for j in 2:Ny
T_v[:, j] .= 0.5 .* (T[:, j-1] .+ T[:, j])
end
T_v[:, 1] .= T_hot # 下壁
T_v[:, Ny+1] .= T_cold # 上壁
@inbounds for j in 2:Ny, i in 1:Nx
v_c = v[i, j]
rho_f = rho_v[i, j] # [kg/m³]
# 対流項 [m/s²]
if i < Nx
u_e = 0.5*(u[i+1, j-1] + u[i+1, j])
v_e = 0.5*(v[i, j] + v[i+1, j])
else
u_e = 0.0; v_e = 0.0
end
if i > 1
u_w = 0.5*(u[i, j-1] + u[i, j])
v_w = 0.5*(v[i-1, j] + v[i, j])
else
u_w = 0.0; v_w = 0.0
end
conv_x = (u_e*v_e - u_w*v_w) / dx
v_n = (j < Ny) ? 0.5*(v[i, j] + v[i, j+1]) : 0.0
v_s = (j > 2) ? 0.5*(v[i, j-1] + v[i, j]) : 0.0
conv_y = (v_n^2 - v_s^2) / dy
# 拡散項 [m/s²]
v_im1 = (i > 1) ? v[i-1, j] : -v[i, j]
v_ip1 = (i < Nx) ? v[i+1, j] : -v[i, j]
diff_x = (v_ip1 - 2*v[i, j] + v_im1) / dx^2
diff_y = (v[i, j+1] - 2*v[i, j] + v[i, j-1]) / dy^2
# 圧力勾配項 [m/s²]
dp_dy = (p[i, j] - p[i, j-1]) / dy # [Pa/m]
# 浮力項：(ρ - ρ_ref) * g / ρ [m/s²]
# 高温（低密度）→ ρ < ρ_ref → 上向き加速（正）
# 低温（高密度）→ ρ > ρ_ref → 下向き加速（負）
rho_v_face = p0 / (R_gas * max(T_v[i, j], 1.0))
buoyancy = -(rho_v_face - rho_ref) * g / rho_f
v_star[i, j] = v_c + dt * (
-(conv_x + conv_y)
+ nu * (diff_x + diff_y)
- dp_dy / rho_f
+ buoyancy
)
end
v_star[:, 1] .= 0.0
v_star[:, Ny+1] .= 0.0
end
# ============================================================
# 発散の計算
# ============================================================
"""速度場の発散を計算 [1/s]"""
function compute_divergence!(div, u, v)
@inbounds for j in 1:Ny, i in 1:Nx
div[i, j] = (u[i+1, j] - u[i, j]) / dx +
(v[i, j+1] - v[i, j]) / dy
end
end
# ============================================================
# 密度変化ソース項
# S_T = α ∇²T / T [1/s]
#
# 導出：
# ∇·u = (1/T) DT/Dt = (1/T)(α∇²T) = α∇²T / T
# ============================================================
"""密度変化ソース項 S_T = α∇²T / T を計算 [1/s]"""
function compute_density_source!(S_T, T)
T_ext = get_T_extended(T)
@inbounds for j in 1:Ny, i in 1:Nx
# ∇²T（中心差分）[K/m²]
lap_T = (T_ext[i+2, j+1] - 2*T_ext[i+1, j+1] + T_ext[i, j+1]) / dx^2 +
(T_ext[i+1, j+2] - 2*T_ext[i+1, j+1] + T_ext[i+1, j ]) / dy^2
# S_T = α ∇²T / T [1/s]
S_T[i, j] = alpha * lap_T / max(T[i, j], 1.0)
end
end
# ============================================================
# 圧力ポアソン方程式（SOR法）
# ∇·(1/ρ ∇p') = (∇·u* - S_T) / dt
# 単位：左辺 [1/(kg/m³) * Pa/m²] = [m/s² / m] = [1/s²]
# 右辺 [(1/s)] / [s] = [1/s²] ✓
# ============================================================
"""圧力修正量 p' を SOR法で求める [Pa]"""
function solve_pressure_poisson!(p_prime, rhs, rho)
fill!(p_prime, 0.0)
for iter in 1:n_pressure_iter
max_res = 0.0
@inbounds for j in 1:Ny, i in 1:Nx
# ノイマン境界：ゴーストセル = 内側セル値
p_e = (i < Nx) ? p_prime[i+1, j] : p_prime[i, j]
p_w = (i > 1) ? p_prime[i-1, j] : p_prime[i, j]
p_n = (j < Ny) ? p_prime[i, j+1] : p_prime[i, j]
p_s = (j > 1) ? p_prime[i, j-1] : p_prime[i, j]
# 面上密度（算術平均）[kg/m³]
rho_e = (i < Nx) ? 0.5*(rho[i, j] + rho[i+1, j]) : rho[i, j]
rho_w = (i > 1) ? 0.5*(rho[i-1, j] + rho[i, j]) : rho[i, j]
rho_n = (j < Ny) ? 0.5*(rho[i, j] + rho[i, j+1]) : rho[i, j]
rho_s = (j > 1) ? 0.5*(rho[i, j-1] + rho[i, j]) : rho[i, j]
# 係数 [m²/kg]
a_e = 1.0 / (rho_e * dx^2)
a_w = 1.0 / (rho_w * dx^2)
a_n = 1.0 / (rho_n * dy^2)
a_s = 1.0 / (rho_s * dy^2)
a_p = a_e + a_w + a_n + a_s
p_new = (a_e*p_e + a_w*p_w + a_n*p_n + a_s*p_s - rhs[i, j]) / a_p
res = abs(p_new - p_prime[i, j])
max_res = max(max_res, res)
p_prime[i, j] = (1.0 - sor_omega)*p_prime[i, j] + sor_omega*p_new
end
# 基準点固定（左下コーナー）
p0_ref = p_prime[1, 1]
p_prime .-= p0_ref
if max_res < pressure_tol
break
end
if iter == n_pressure_iter
println("res= ", max_res)
end
end
end
# ============================================================
# 温度方��式（陽的，風上差分＋中心差分）
# ∂T/∂t + u·∇T = α ∇²T
# ============================================================
"""温度方程式を1ステップ進める [K]"""
function advance_temperature!(T_new, u, v, T)
T_ext = get_T_extended(T)
@inbounds for j in 1:Ny, i in 1:Nx
u_c = 0.5*(u[i, j] + u[i+1, j]) # [m/s]
v_c = 0.5*(v[i, j] + v[i, j+1]) # [m/s]
# 対流項（1次風上差分）[K/s]
dTdx_p = (T_ext[i+2, j+1] - T_ext[i+1, j+1]) / dx
dTdx_m = (T_ext[i+1, j+1] - T_ext[i, j+1]) / dx
dTdy_p = (T_ext[i+1, j+2] - T_ext[i+1, j+1]) / dy
dTdy_m = (T_ext[i+1, j+1] - T_ext[i+1, j ]) / dy
conv = u_c * (u_c >= 0 ? dTdx_m : dTdx_p) +
v_c * (v_c >= 0 ? dTdy_m : dTdy_p)
# 拡散項（中心差分）[K/s]
diff = alpha * (
(T_ext[i+2, j+1] - 2*T_ext[i+1, j+1] + T_ext[i, j+1]) / dx^2 +
(T_ext[i+1, j+2] - 2*T_ext[i+1, j+1] + T_ext[i+1, j ]) / dy^2
)
T_new[i, j] = T[i, j] + dt * (-conv + diff)
end
end
# ============================================================
# PISOアルゴリズム（1ステップ）
# ============================================================
"""PISOアルゴリズムで1タイムステップ進める"""
function piso_step!(u, v, p, T, rho,
u_star, v_star, p_prime, div, S_T, T_new)
# ステップ0：密度の更新 [kg/m³]
compute_density!(rho, T)
rho_u = interpolate_rho_u(rho)
rho_v = interpolate_rho_v(rho)
# ステップ1：中間速度（u*, v*）の計算 [m/s]
compute_u_star!(u_star, u, v, p, rho_u)
compute_v_star!(v_star, u, v, p, T, rho_v)
apply_velocity_bc!(u_star, v_star)
# 密度変化ソース項 [1/s]
compute_density_source!(S_T, T)
# ステップ2：圧力修正ループ（PISO）
u_c = copy(u_star)
v_c = copy(v_star)
for _ in 1:n_correctors
compute_divergence!(div, u_c, v_c)
# 右辺 [1/s²]
rhs = (div .- S_T) ./ dt
solve_pressure_poisson!(p_prime, rhs, rho)
p .+= p_prime
# 速度修正：u -= dt/ρ * ∇p' [m/s]
@inbounds for j in 1:Ny, i in 2:Nx
u_c[i, j] -= (dt / rho_u[i, j]) * (p_prime[i, j] - p_prime[i-1, j]) / dx
end
@inbounds for j in 2:Ny, i in 1:Nx
v_c[i, j] -= (dt / rho_v[i, j]) * (p_prime[i, j] - p_prime[i, j-1]) / dy
end
apply_velocity_bc!(u_c, v_c)
end
# 動的圧力の基準値固定（平均ゼロ）
p .-= mean(p)
u .= u_c
v .= v_c
# ステップ3：温度方程式を進める
advance_temperature!(T_new, u, v, T)
T .= T_new
end
# ============================================================
# ヌッセルト数の計算
# Nu = 1 + <v T> / (α ΔT / Ly)
# ============================================================
"""体積平均ヌッセルト数を計算"""
function compute_nusselt(T, v)
v_c = 0.5 .* (v[:, 1:Ny] .+ v[:, 2:Ny+1]) # [m/s]
T_nd = (T .- T_cold) ./ (T_hot - T_cold) # 無次元温度（表示用）
Nu = 1.0 + mean(v_c .* T_nd) * Ly / alpha
return Nu
end
# ============================================================
# 結果の出力（CSV）
# ============================================================
"""温度場・速度場・密度場をCSVに出力"""
function save_fields(u, v, p, T, rho, step)
filename = @sprintf("rb_lowmach_dim_step%06d.csv", step)
u_c = 0.5 .* (u[1:Nx, :] .+ u[2:Nx+1, :])
v_c = 0.5 .* (v[:, 1:Ny] .+ v[:, 2:Ny+1])
open(filename, "w") do io
println(io, "x[m],y[m],u[m/s],v[m/s],p[Pa],T[K],rho[kg/m3]")
for j in 1:Ny, i in 1:Nx
@printf(io, "%.6f,%.6f,%.6e,%.6e,%.6e,%.4f,%.6f\n",
xc[i], yc[j], u_c[i,j], v_c[i,j], p[i,j], T[i,j], rho[i,j])
end
end
end
"""ヌッセルト数履歴をCSVに出力"""
function save_nusselt(t_history, Nu_history)
open("nusselt_history.csv", "w") do io
println(io, "t[s],Nu")
for (t, Nu) in zip(t_history, Nu_history)
@printf(io, "%.4f,%.6f\n", t, Nu)
end
end
end
# ============================================================
# メインループ
# ============================================================
function main()
println("=" ^ 65)
println("2D Rayleigh-Bénard Convection (Low-Mach, Dimensional)")
@printf(" Ra = %.3e, Pr = %.3f\n", Ra, Pr)
@printf(" Nx = %d, Ny = %d, Lx = %.1f m, Ly = %.1f m\n", Nx, Ny, Lx, Ly)
@printf(" dx = %.4f m, dy = %.4f m\n", dx, dy)
@printf(" dt = %.2e s, t_end = %.1f s\n", dt, t_end)
println("-" ^ 65)
@printf(" T_hot = %.1f K, T_cold = %.1f K, T_init = %.1f K\n", T_hot, T_cold, T_init)
@printf(" T_ref = %.1f K, p0 = %.1f Pa\n", T_ref, p0)
@printf(" ρ_ref = %.4f kg/m³\n", rho_ref)
@printf(" ν = %.3e m²/s\n", nu)
@printf(" α = %.3e m²/s\n", alpha)
@printf(" λ = %.4f W/(m·K)\n", lambda)
println("=" ^ 65)
# 安定性チェック
u_est = sqrt(g * (1.0/T_ref) * (T_hot - T_cold) * Ly) # 自由対流速度スケール
cfl = u_est * dt / min(dx, dy)
diff_num = max(nu, alpha) * dt / min(dx, dy)^2
@printf(" u_scale = %.4f m/s\n", u_est)
@printf(" CFL(est) = %.3f, Diffusion number = %.3f\n", cfl, diff_num)
if cfl > 0.5 || diff_num > 0.5
println(" WARNING: Time step may be too large for stability!")
end
println()
# 変数の初期化
u = zeros(Nx+1, Ny) # [m/s]
v = zeros(Nx, Ny+1) # [m/s]
p = zeros(Nx, Ny) # 動的圧力 [Pa]
T = fill(T_init, Nx, Ny) # [K]
rho = zeros(Nx, Ny) # [kg/m³]
# 初期温度場：T_init = T_cold = 300 K + 微小擾乱
#for j in 1:Ny, i in 1:Nx
# T[i, j] = T_init + 0.1 * randn() # [K]，擾乱スケール 0.1 K
#end
#clamp!(T, T_cold - 5.0, T_hot + 5.0)
compute_density!(rho, T)
# 作業バッファ
u_star = zeros(Nx+1, Ny)
v_star = zeros(Nx, Ny+1)
p_prime = zeros(Nx, Ny)
div = zeros(Nx, Ny)
S_T = zeros(Nx, Ny)
T_new = zeros(Nx, Ny)
t = 0.0
step = 0
Nu_history = Float64[]
t_history = Float64[]
# ウォームアップ（JIT）
#println(" JIT warmup...")
#for _ in 1:5
# piso_step!(u, v, p, T, rho, u_star, v_star, p_prime, div, S_T, T_new)
#end
# リセット
fill!(u, 0.0); fill!(v, 0.0); fill!(p, 0.0)
#for j in 1:Ny, i in 1:Nx
# T[i, j] = T_init + 0.1 * randn()
#end
#clamp!(T, T_cold - 5.0, T_hot + 5.0)
compute_density!(rho, T)
t = 0.0; step = 0
println(" Starting main loop...")
t_cpu_start = time()
while t < t_end
piso_step!(u, v, p, T, rho, u_star, v_star, p_prime, div, S_T, T_new)
t += dt
step += 1
if step % save_interval == 0
Nu = compute_nusselt(T, v)
rho_min = minimum(rho)
rho_max = maximum(rho)
T_min = minimum(T)
T_max = maximum(T)
u_max = maximum(abs.(u))
v_max = maximum(abs.(v))
elapsed = time() - t_cpu_start
push!(Nu_history, Nu)
push!(t_history, t)
@printf(" step=%6d t=%8.2f s Nu=%7.4f |u|max=%7.4f m/s |v|max=%7.4f m/s\n",
step, t, Nu, u_max, v_max)
@printf(" T=[%6.2f,%6.2f] K ρ=[%.4f,%.4f] kg/m³ CPU=%.1f s\n",
T_min, T_max, rho_min, rho_max, elapsed)
save_fields(u, v, p, T, rho, step)
end
end
save_nusselt(t_history, Nu_history)
println("\nSimulation finished.")
@printf(" Total CPU time : %.2f s\n", time() - t_cpu_start)
@printf(" Final Nu : %.4f\n", Nu_history[end])
println(" Saved: rb_lowmach_dim_step*.csv, nusselt_history.csv")
end
main()
Juliaのソースコードから説明の文章もTeXで作ってもらいました。こっちはWebのClaudeでやりましたが，これもすんなりコンパイルできるものが生成できました。
[Read More]
