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でやりましたが，これもすんなりコンパイルできるものが生成できました。
TeXのソース （ここをクリックすると詳細が出ます）
%% rayleigh_benard_lowmach_description.tex
%% Engine: LuaLaTeX / Class: ltjsarticle
\documentclass[a4paper,10pt]{ltjsarticle}
%% ---------------------------------------------------------------
%% パッケージ
%% ---------------------------------------------------------------
\usepackage{luatexja}
\usepackage{luatexja-fontspec}
\usepackage{amsmath,amssymb,mathtools}
\usepackage{booktabs}
\usepackage{array}
\usepackage{geometry}
\usepackage{hyperref}
\usepackage{xcolor}
\usepackage{enumitem}
\usepackage{listings}
\usepackage{fancyhdr}
\geometry{top=25mm, bottom=25mm, left=25mm, right=25mm}
%% ---------------------------------------------------------------
%% ヘッダ・フッタ
%% ---------------------------------------------------------------
\pagestyle{fancy}
\fancyhf{}
\fancyhead[L]{\small 2次元レイリー・ベナール対流解析コード（低マッハ数近似・次元あり）}
\fancyhead[R]{\small \thepage}
\renewcommand{\headrulewidth}{0.4pt}
%% ---------------------------------------------------------------
%% listings 設定（Julia）
%% ---------------------------------------------------------------
\lstdefinelanguage{Julia}{
keywords={function,end,for,if,else,elseif,return,const,using,while,in,do,push!,fill!,copy,zeros,ones,println,break},
sensitive=true,
comment=[l]{\#},
string=[b]",
morestring=[b]',
}
\lstset{
language=Julia,
basicstyle=\ttfamily\small,
keywordstyle=\color{blue}\bfseries,
commentstyle=\color{gray}\itshape,
stringstyle=\color{teal},
numbers=left,
numberstyle=\tiny\color{gray},
frame=single,
breaklines=true,
tabsize=4,
}
%% ---------------------------------------------------------------
%% 数式マクロ
%% ---------------------------------------------------------------
\newcommand{\bnabla}{\boldsymbol{\nabla}}
\newcommand{\bu}{\boldsymbol{u}}
\newcommand{\bv}{\boldsymbol{v}}
\newcommand{\DP}[2]{\dfrac{\partial #1}{\partial #2}}
\newcommand{\DDt}[1]{\dfrac{D #1}{D t}}
\newcommand{\rhoref}{\rho_{\mathrm{ref}}}
\newcommand{\Tref}{T_{\mathrm{ref}}}
%% ---------------------------------------------------------------
\begin{document}
\title{%
\textbf{2次元レイリー・ベナール対流解析コード}\\
\large 低マッハ数近似・次元あり（Julia実装）\\[1ex]
\normalsize \texttt{rayleigh\_benard\_lowmach\_dimensional.jl} 処理内容の説明
}
\author{}
\date{}
\maketitle
\thispagestyle{fancy}
\tableofcontents
\newpage
%% ===============================================================
\section{概要}
%% ===============================================================
本コードは，密度差による浮力が駆動力となる自然対流現象，
\textbf{レイリー・ベナール対流}（Rayleigh--B\'{e}nard Convection）を
2次元で数値シミュレーションする Julia プログラムである．
作動流体として空気（理想気体）を用い，
\textbf{低マッハ数近似}（Low-Mach Number Approximation）のもとで
\textbf{SI単位系の次元あり}物理量を扱う．
数値アルゴリズムの主な特徴を以下に示す．
\begin{itemize}
\item \textbf{低マッハ数近似}：熱力学的圧力 $p_0$（空間一様・時間変動なし）と
動的圧力 $p'$（速度・密度変動に寄与）を分離する．
密度は状態方程式 $\rho = p_0/(R_\mathrm{gas}\,T)$ で決定される．
\item \textbf{PISO法}（Pressure Implicit with Splitting of Operators）：
圧力--速度連成を反復的に解くための分離型アルゴリズム．
\item \textbf{スタガード格子}（Staggered Grid）：
$u$，$v$，$T$，$p$ をそれぞれ異なる格子点に配置し，
速度--圧力の空間的なデカップリングを防ぐ．
\item 境界条件：全壁すべりなし（no-slip），
下壁高温（$T_\mathrm{hot}$），上壁低温（$T_\mathrm{cold}$），
側壁断熱（Neumann）．
\end{itemize}
%% ===============================================================
\section{支配方程式}
%% ===============================================================
\subsection{低マッハ数近似の連続方程式}
低マッハ数近似では音波を除去するため，
熱力学的圧力 $p_0$ を空間一様定数として扱う．
理想気体の状態方程式は
\begin{equation}
\rho = \frac{p_0}{R_\mathrm{gas}\,T}
\label{eq:eos}
\end{equation}
となる．連続方程式（質量保存）は
\begin{equation}
\bnabla \cdot \bu = S_T \equiv \frac{\alpha\,\nabla^2 T}{T}
\label{eq:cont}
\end{equation}
と非ゼロの発散源項 $S_T$ をもつ形になる．
これは温度変化に伴う密度変化を表す．
\subsection{運動量方程式}
$x$ 方向（水平），$y$ 方向（鉛直上向き）の運動量方程式を以下に示す．
\begin{align}
\rho\,\DDt{u} &= -\DP{p'}{x} + \mu\,\nabla^2 u
\label{eq:mom_x} \\[6pt]
\rho\,\DDt{v} &= -\DP{p'}{y} + \mu\,\nabla^2 v
+ (\rho - \rhoref)\,g
\label{eq:mom_y}
\end{align}
ここで $\mu$ は動粘性率，$g$ は重力加速度，$\rhoref$ は参照密度である．
式~\eqref{eq:mom_y} 右辺第3項が\textbf{浮力項}（Boussinesq 型）であり，
高温低密度域では上向き，低温高密度域では下向きの加速を生じる．
\subsection{エネルギー方程式（温度方程式）}
\begin{equation}
\DP{T}{t} + \bu\cdot\bnabla T = \alpha\,\nabla^2 T
\label{eq:energy}
\end{equation}
ここで $\alpha$ は熱拡散率（定数）である．
\subsection{無次元パラメータ}
\begin{align}
\mathrm{Ra} &= \frac{g\,\beta\,(T_\mathrm{hot}-T_\mathrm{cold})\,L_y^3}{\nu\,\alpha}
\label{eq:Ra} \\[4pt]
\mathrm{Pr} &= \frac{\nu}{\alpha}
\label{eq:Pr}
\end{align}
ここで $\beta = 1/\Tref$（理想気体の体積膨張係数），
$\nu = \mu/\rhoref$ は動粘性係数，$L_y$ は鉛直方向の領域高さである．
%% ===============================================================
\section{物理定数と計算条件}
%% ===============================================================
\begin{table}[htbp]
\centering
\caption{物理定数・流体物性（空気，参照温度 $\Tref = 375\,\mathrm{K}$）}
\label{tab:constants}
\begin{tabular}{llll}
\toprule
パラメータ & 記号 & 値 & 単位 \\
\midrule
重力加速度 & $g$ & $9.81$ & $\mathrm{m/s^2}$ \\
気体定数（空気） & $R_\mathrm{gas}$ & $287.0$ & $\mathrm{J/(kg\cdot K)}$ \\
定圧比熱 & $c_p$ & $1005.0$ & $\mathrm{J/(kg\cdot K)}$ \\
下壁温度 & $T_\mathrm{hot}$ & $450.0$ & $\mathrm{K}$ \\
上壁温度 & $T_\mathrm{cold}$ & $300.0$ & $\mathrm{K}$ \\
熱力学的圧力 & $p_0$ & $101325.0$ & $\mathrm{Pa}$ \\
参照密度 & $\rhoref$ & $p_0/(R_\mathrm{gas}\,\Tref)\approx 1.084$ & $\mathrm{kg/m^3}$ \\
動粘性係数 & $\nu$ & $1.0\times10^{-2}$ & $\mathrm{m^2/s}$ \\
熱拡散率 & $\alpha$ & $1.0\times10^{-2}$ & $\mathrm{m^2/s}$ \\
熱伝導率 & $\lambda$ & $0.0281$ & $\mathrm{W/(m\cdot K)}$ \\
\bottomrule
\end{tabular}
\end{table}
\begin{table}[htbp]
\centering
\caption{計算領域・格子・時間ステップ}
\label{tab:grid}
\begin{tabular}{llll}
\toprule
パラメータ & 記号 & 値 & 単位 \\
\midrule
$x$方向領域長さ & $L_x$ & $2.0$ & $\mathrm{m}$ \\
$y$方向領域長さ & $L_y$ & $1.0$ & $\mathrm{m}$ \\
$x$方向セル数 & $N_x$ & $64$ & --- \\
$y$方向セル数 & $N_y$ & $32$ & --- \\
時間刻み & $\Delta t$ & $5\times10^{-4}$ & $\mathrm{s}$ \\
終了時間 & $t_\mathrm{end}$ & $10.0$ & $\mathrm{s}$ \\
PISO修正回数 & $n_\mathrm{corr}$ & $2$ & --- \\
圧力 SOR 反復数 & --- & $1500$ & --- \\
SOR 緩和係数 & $\omega$ & $1.6$ & --- \\
\bottomrule
\end{tabular}
\end{table}
%% ===============================================================
\section{スタガード格子の配置}
%% ===============================================================
本コードでは，圧力--速度のデカップリングを防ぐためにスタガード格子を採用している
（図\ref{fig:staggered} 参照）．
\begin{figure}[htbp]
\centering
\setlength{\unitlength}{1.0mm}
\begin{picture}(90,70)
%% セル境界
\put(10,10){\line(1,0){70}}
\put(10,40){\line(1,0){70}}
\put(10,70){\line(1,0){70}}
\put(10,10){\line(0,1){60}}
\put(45,10){\line(0,1){60}}
\put(80,10){\line(0,1){60}}
%% u 速度（左右面中央）
\put(10,25){\circle*{2}} \put(4,23){\small$u_{i,j}$}
\put(45,25){\circle*{2}} \put(39,23){\small$u_{i+1,j}$}
\put(80,25){\circle*{2}}
\put(10,55){\circle*{2}}
\put(45,55){\circle*{2}}
\put(80,55){\circle*{2}}
%% v 速度（上下面中央）
\put(27,10){\circle{3}} \put(23,5){\small$v_{i,j}$}
\put(62,10){\circle{3}}
\put(27,40){\circle{3}} \put(23,35){\small$v_{i,j+1}$}
\put(62,40){\circle{3}}
\put(27,70){\circle{3}}
\put(62,70){\circle{3}}
%% p, T（セル中央）
\put(27,25){\circle{4}} \put(20,22){\small$p,T$}
\put(62,25){\circle{4}}
\put(27,55){\circle{4}}
\put(62,55){\circle{4}}
%% ラベル
\put(23,28){\small$(i,j)$}
\put(56,28){\small$(i{+}1,j)$}
\put(50,0){\small $x$}
\put(0,50){\small $y$}
\end{picture}
\caption{スタガード格子の変数配置．
$\bullet$：$u$（$x$方向速度，セル左右面中央），
$\circ$（中空）：$v$（$y$方向速度，セル上下面中央），
$\circ$（太）：$p$，$T$，$\rho$（セル中央）．}
\label{fig:staggered}
\end{figure}
各変数の配列サイズは以下の通りである．
\begin{itemize}
\item $u$：$(N_x+1) \times N_y$ （$x$面，$i=1,\ldots,N_x+1$）
\item $v$：$N_x \times (N_y+1)$ （$y$面，$j=1,\ldots,N_y+1$）
\item $p$，$T$，$\rho$：$N_x \times N_y$ （セル中央）
\end{itemize}
%% ===============================================================
\section{境界条件}
%% ===============================================================
\subsection{速度境界条件}
全壁面でのすべりなし・不透過条件を課す：
\begin{equation}
u\big|_{i=1,\,N_x+1} = 0, \qquad
v\big|_{j=1,\,N_y+1} = 0.
\end{equation}
\subsection{温度境界条件}
温度のゴーストセル法により壁境界条件を実装する：
\begin{align}
&\text{下壁（Dirichlet，高温）：} && T_\mathrm{ghost} = 2\,T_\mathrm{hot} - T_{i,1} \\
&\text{上壁（Dirichlet，低温）：} && T_\mathrm{ghost} = 2\,T_\mathrm{cold} - T_{i,N_y} \\
&\text{左右壁（Neumann，断熱）：} && T_\mathrm{ghost} = T_\mathrm{boundary}
\end{align}
%% ===============================================================
\section{PISOアルゴリズムの流れ}
%% ===============================================================
1タイムステップの処理は以下の手順（\texttt{piso\_step!} 関数）で実行される．
\subsection*{ステップ 0：密度の更新}
理想気体の状態方程式~\eqref{eq:eos} を用いて全セルの密度を計算する．
続いて，$u$面・$v$面上への線形補間によりそれぞれの面中央密度
$\rho_u$，$\rho_v$ を求める：
\begin{equation}
\rho_{u,i,j} = \frac{\rho_{i-1,j} + \rho_{i,j}}{2}, \qquad
\rho_{v,i,j} = \frac{\rho_{i,j-1} + \rho_{i,j}}{2}.
\end{equation}
\subsection*{ステップ 1：中間速度の計算}
陽的 Euler 法により中間速度（$u^*$，$v^*$）を計算する：
\begin{align}
u^* &= u^n + \Delta t\left[
-(\text{対流項}_x)
+ \nu\,\nabla^2 u^n
- \frac{1}{\rho_u}\DP{p^n}{x}
\right] \\[6pt]
v^* &= v^n + \Delta t\left[
-(\text{対流項}_y)
+ \nu\,\nabla^2 v^n
- \frac{1}{\rho_v}\DP{p^n}{y}
+ \frac{(\rhoref - \rho_v)\,g}{\rho_v}
\right]
\end{align}
対流項は中心差分，拡散項は2階中心差分で離散化される．
$u$ 方程式は内壁スタガード点 $i=2,\ldots,N_x$，
$v$ 方程式は $j=2,\ldots,N_y$ で計算し，
壁面値はゼロ（no-slip）を直接設定する．
\subsection*{ステップ 2：圧力修正ループ（PISO）}
PISO 法では圧力修正を $n_\mathrm{corr}$ 回繰り返す（本コードでは2回）．
\paragraph{発散の計算}
\begin{equation}
(\bnabla\cdot\bu^*)_{i,j}
= \frac{u^*_{i+1,j} - u^*_{i,j}}{\Delta x}
+ \frac{v^*_{i,j+1} - v^*_{i,j}}{\Delta y}
\end{equation}
\paragraph{密度変化ソース項}
\begin{equation}
S_{T,i,j} = \frac{\alpha\,\nabla^2 T}{T}\bigg|_{i,j}
\label{eq:ST}
\end{equation}
これは連続方程式~\eqref{eq:cont} の非ゼロ発散源項に対応する．
\paragraph{圧力ポアソン方程式}
\begin{equation}
\bnabla\cdot\!\left(\frac{1}{\rho}\,\bnabla p'\right)
= \frac{\bnabla\cdot\bu^* - S_T}{\Delta t}
\label{eq:poisson}
\end{equation}
式~\eqref{eq:poisson} を，SOR 法（逐次過緩和法）により反復求解する：
\begin{equation}
p'^{(k+1)}_{i,j}
= (1-\omega)\,p'^{(k)}_{i,j}
+ \omega\,\frac{a_e\,p'_e + a_w\,p'_w + a_n\,p'_n + a_s\,p'_s - b_{i,j}}{a_e+a_w+a_n+a_s}
\end{equation}
ここで係数は
$a_e = 1/(\rho_e\,\Delta x^2)$，
$a_w = 1/(\rho_w\,\Delta x^2)$，
$a_n = 1/(\rho_n\,\Delta y^2)$，
$a_s = 1/(\rho_s\,\Delta y^2)$
であり，$\omega = 1.6$ は SOR 緩和係数である．
全壁面での圧力修正量境界条件は Neumann 条件（$\partial p'/\partial n=0$）とし，
左下コーナーを基準点（$p'=0$）として一意性を確保する．
\paragraph{速度・圧力の修正}
\begin{align}
u_{i,j} &\leftarrow u^*_{i,j}
- \frac{\Delta t}{\rho_{u,i,j}}\,\frac{p'_{i,j} - p'_{i-1,j}}{\Delta x} \\[4pt]
v_{i,j} &\leftarrow v^*_{i,j}
- \frac{\Delta t}{\rho_{v,i,j}}\,\frac{p'_{i,j} - p'_{i,j-1}}{\Delta y} \\[4pt]
p &\leftarrow p + p'
\end{align}
動的圧力 $p$ は空間平均ゼロになるよう毎ステップ基準化する．
\subsection*{ステップ 3：温度方程式の前進}
式~\eqref{eq:energy} を陽的1次精度風上差分で離散化する：
\begin{equation}
T^{n+1}_{i,j} = T^n_{i,j}
+ \Delta t\left[
-u_c\,\delta_x T - v_c\,\delta_y T
+ \alpha\,\nabla^2 T
\right]^n_{i,j}
\end{equation}
対流項は速度の符号に応じた風上差分（1次精度），
拡散項は中心差分（2次精度）を用いる．
$u_c$，$v_c$ はセル中央速度（スタガード速度の線形補間値）である．
\subsection*{PISOアルゴリズムのフロー}
\begin{enumerate}[label=\textbf{Step~\arabic*.}, leftmargin=*]
\item 状態方程式 $\rho = p_0/(R_\mathrm{gas}T)$ で密度更新，面上密度を補間
\item 中間速度 $(u^*, v^*)$ を陽的に計算
\item $n_\mathrm{corr}$ 回のPISO修正ループ：
\begin{enumerate}[label=\alph*.]
\item 発散 $\bnabla\cdot\bu^*$ を計算
\item 密度ソース $S_T$ を計算
\item 圧力ポアソン方程式をSORで解き $p'$ を得る
\item $p \leftarrow p + p'$，速度を $p'$ で修正
\end{enumerate}
\item 動的圧力を平均ゼロに基準化
\item 温度方程式を1ステップ前進
\end{enumerate}
%% ===============================================================
\section{診断量：ヌッセルト数}
%% ===============================================================
対流伝熱の強さを表す体積平均ヌッセルト数を以下の式で計算する：
\begin{equation}
\mathrm{Nu} = 1 + \frac{\langle v'\,\theta\rangle}{\alpha\,\Delta T / L_y}
\label{eq:Nu}
\end{equation}
ここで $\theta = (T - T_\mathrm{cold})/(T_\mathrm{hot}-T_\mathrm{cold})$ は無次元温度，
$\langle\cdot\rangle$ は体積平均，
$\Delta T = T_\mathrm{hot} - T_\mathrm{cold}$ である．
$\mathrm{Nu}=1$ は純粋な熱伝導，$\mathrm{Nu}>1$ が対流伝熱の寄与を示す．
%% ===============================================================
\section{出力}
%% ===============================================================
\subsection{物理場のCSV出力（\texttt{save\_fields}）}
$n_\mathrm{save} = 1000$ ステップごとに物理場スナップショットをCSVファイルに保存する．
\begin{center}
\texttt{rb\_lowmach\_dim\_step\{XXXXXX\}.csv}
\end{center}
出力される列は以下の7項目である：
\begin{center}
\begin{tabular}{ll}
\toprule
列名 & 内容 \\
\midrule
\texttt{x[m]} & $x$座標（セル中央）[m] \\
\texttt{y[m]} & $y$座標（セル中央）[m] \\
\texttt{u[m/s]} & $x$方向速度（セル中央補間値）[m/s] \\
\texttt{v[m/s]} & $y$方向速度（セル中央補間値）[m/s] \\
\texttt{p[Pa]} & 動的圧力 [Pa] \\
\texttt{T[K]} & 温度 [K] \\
\texttt{rho[kg/m3]} & 密度 [kg/m$^3$] \\
\bottomrule
\end{tabular}
\end{center}
\subsection{ヌッセルト数履歴（\texttt{save\_nusselt}）}
計算終了後に \texttt{nusselt\_history.csv} としてヌッセルト数の時系列を保存する．
\subsection{コンソール出力}
各保存ステップでは以下の情報を標準出力に表示する：
\begin{itemize}
\item ステップ数・シミュレーション時刻
\item ヌッセルト数 $\mathrm{Nu}$
\item 最大速度 $|u|_\mathrm{max}$，$|v|_\mathrm{max}$
\item 温度の最小・最大値
\item 密度の最小・最大値
\item 経過CPU時間
\end{itemize}
%% ===============================================================
\section{安定性条件のチェック}
%% ===============================================================
計算開始前に，推定自由対流速度スケール
\begin{equation}
u_\mathrm{scale} = \sqrt{g\,\beta\,\Delta T\,L_y}
\end{equation}
をもとに，CFL数および拡散数を評価して警告を表示する：
\begin{align}
C &= \frac{u_\mathrm{scale}\,\Delta t}{\min(\Delta x,\Delta y)} \leq 0.5 \\[4pt]
D &= \frac{\max(\nu,\alpha)\,\Delta t}{\min(\Delta x,\Delta y)^2} \leq 0.5
\end{align}
%% ===============================================================
\section{主な関数の一覧}
%% ===============================================================
\begin{table}[htbp]
\centering
\caption{主な関数と役割}
\label{tab:functions}
\begin{tabular}{>{\ttfamily}lp{9cm}}
\toprule
\normalfont 関数名 & 役割 \\
\midrule
compute\_density! & 理想気体状態方程式 $\rho=p_0/(R_\mathrm{gas}T)$ で密度更新 \\
interpolate\_rho\_u & $u$面上への密度の線形補間 \\
interpolate\_rho\_v & $v$面上への密度の線形補間 \\
apply\_velocity\_bc! & すべりなし速度境界条件の適用 \\
get\_T\_extended & 温度ゴーストセルの生成（壁温度・断熱境界） \\
compute\_u\_star! & 中間速度 $u^*$ の計算（対流・拡散・圧力勾配項） \\
compute\_v\_star! & 中間速度 $v^*$ の計算（浮力項を含む） \\
compute\_divergence! & 速度発散 $\bnabla\cdot\bu$ の計算 \\
compute\_density\_source! & 密度変化ソース項 $S_T = \alpha\nabla^2T/T$ の計算 \\
solve\_pressure\_poisson! & SOR法による圧力ポアソン方程式の求解 \\
advance\_temperature! & 温度方程式の陽的前進（風上＋中心差分） \\
piso\_step! & PISOアルゴリズム1ステップの統合処理 \\
compute\_nusselt & 体積平均ヌッセルト数の計算 \\
save\_fields & 物理場スナップショットのCSV出力 \\
save\_nusselt & ヌッセルト数時系列のCSV出力 \\
main & 初期化・メインループ・結果出力 \\
\bottomrule
\end{tabular}
\end{table}
%% ===============================================================
\section{実行方法}
%% ===============================================================
Julia の REPL または以下のコマンドで実行する：
\begin{lstlisting}[language=bash, numbers=none]
julia rayleigh_benard_lowmach_dimensional.jl
\end{lstlisting}
\noindent 依存パッケージは標準ライブラリのみ（\texttt{Printf}，\texttt{LinearAlgebra}，\texttt{Statistics}）であり，
追加インストールは不要である．
%% ===============================================================
\appendix
\section{密度変化ソース項 \texorpdfstring{$S_T$}{S\_T} の導出}
\label{app:ST_derivation}
%% ===============================================================
\subsection{出発点：一般の圧縮性連続方程式}
質量保存則（圧縮性連続方程式）は
\begin{equation}
\DP{\rho}{t} + \bnabla\cdot(\rho\,\bu) = 0
\label{eq:app_cont_comp}
\end{equation}
である．これを展開すると
\begin{equation}
\DP{\rho}{t} + \bu\cdot\bnabla\rho + \rho\,\bnabla\cdot\bu = 0
\end{equation}
となり，実質微分 $D/Dt \equiv \partial/\partial t + \bu\cdot\bnabla$ を用いると
\begin{equation}
\frac{1}{\rho}\,\DDt{\rho} + \bnabla\cdot\bu = 0.
\label{eq:app_cont_lagrange}
\end{equation}
\subsection{低マッハ数近似と状態方程式}
低マッハ数近似では，熱力学的圧力 $p_0$ を空間一様・時間定数とみなす．
理想気体の状態方程式は
\begin{equation}
\rho = \frac{p_0}{R_\mathrm{gas}\,T}
\label{eq:app_eos}
\end{equation}
である．$p_0$ は定数であるから，$\rho$ は $T$ のみの関数となる：
\begin{equation}
\rho = \rho(T).
\end{equation}
\subsection{密度の実質微分を温度で表す}
式~\eqref{eq:app_eos} を時間微分する（流体粒子に追随した実質微分）：
\begin{equation}
\DDt{\rho} = \frac{D}{Dt}\!\left(\frac{p_0}{R_\mathrm{gas}\,T}\right)
= -\frac{p_0}{R_\mathrm{gas}\,T^2}\,\DDt{T}.
\label{eq:app_Drho}
\end{equation}
状態方程式~\eqref{eq:app_eos} を用いて $p_0/(R_\mathrm{gas}\,T) = \rho$ と置き換えると
\begin{equation}
\DDt{\rho} = -\frac{\rho}{T}\,\DDt{T}.
\label{eq:app_Drho2}
\end{equation}
これを式~\eqref{eq:app_cont_lagrange} に代入すると
\begin{equation}
-\frac{1}{T}\,\DDt{T} + \bnabla\cdot\bu = 0,
\end{equation}
すなわち
\begin{equation}
\bnabla\cdot\bu = \frac{1}{T}\,\DDt{T}.
\label{eq:app_div_DT}
\end{equation}
\subsection{温度方程式の代入}
エネルギー方程式（温度方程式）は，粘性散逸・圧縮仕事を無視した低マッハ数近似のもとで
\begin{equation}
\DDt{T} = \alpha\,\nabla^2 T
\label{eq:app_energy}
\end{equation}
と書ける（$\alpha$ は熱拡散率）．
式~\eqref{eq:app_div_DT} に式~\eqref{eq:app_energy} を代入すると，
\begin{equation}
\boxed{
\bnabla\cdot\bu
= \frac{1}{T}\,\alpha\,\nabla^2 T
= \frac{\alpha\,\nabla^2 T}{T}
\equiv S_T
}
\label{eq:app_ST_result}
\end{equation}
が得られる．これが密度変化ソース項 $S_T$ の定義式である．
\subsection{物理的意味}
式~\eqref{eq:app_ST_result} は，\textbf{温度の拡散（ラプラシアン）が正の領域では流体が膨張し，
負の領域では収縮する}ことを意味する．
具体的には以下のように解釈できる：
\begin{itemize}
\item $\nabla^2 T > 0$（周囲より低温な領域，熱が流入）：
温度上昇 $\Rightarrow$ 密度低下 $\Rightarrow$ $\bnabla\cdot\bu > 0$（膨張）
\item $\nabla^2 T < 0$（周囲より高温な領域，熱が流出）：
温度低下 $\Rightarrow$ 密度上昇 $\Rightarrow$ $\bnabla\cdot\bu < 0$（収縮）
\end{itemize}
非圧縮性近似（$\rho = \mathrm{const}$）では $S_T \equiv 0$ となり，
$\bnabla\cdot\bu = 0$ の通常の非圧縮条件に戻る．
低マッハ数近似はこの $S_T \neq 0$ を許容することで
\textbf{密度変化を伴う熱対流}を音波を含まずに正確に記述する．
\subsection{圧力ポアソン方程式への組み込み}
速度修正後に $\bnabla\cdot\bu = S_T$ が成立するよう，
圧力修正量 $p'$ のポアソン方程式は
\begin{equation}
\bnabla\cdot\!\left(\frac{1}{\rho}\,\bnabla p'\right)
= \frac{\bnabla\cdot\bu^* - S_T}{\Delta t}
\end{equation}
と修正される（本文式~\eqref{eq:poisson} 参照）．
右辺の $-S_T$ が低マッハ数近似特有の補正項であり，
これを省略すると密度変化の効果が無視されて質量保存が破れる．
\end{document}
モデルとか手法の名前が分かってれば，生成AIでコードを作成することが可能な時代になりました。ただ，その手法のことを分かってないと，合ってるかどうかも全然わからないので，そこの勉強は必要だと思います。
生成AIでコードを作るのはいいと思います。ただ，数値解析のコードは動けばいいってもんじゃなくて，その結果が物理的に正しいかをチェックする必要があります。別件でWebアプリも作りましたが，そっちはただ動けばいいので，特に問題なかったですが，数値解析は正しいか考える必要があります。例えば，境界条件を設定している場所はあってるかとか，境界条件の内容が正しいかとか。
プログラムは人間が書く必要ない気もしてましたが，最近みたネットの情報で，「人間はプログラムのできあがりを管理するプロマネみたいなもん」ってのがありました。（かなり個人的な解釈が入ってます）
元 https://x.com/kaityo256/status/2031701255006171257
プロマネが内容分かってないわけにはいかないので，結局書くのはAIだけど，人間側には内容をちゃんと分かってる能力が必要です。なので，プログラムの勉強はちゃんとする必要があるってのがいまのところの結論です。
プログラムの勉強
これからの時代の賢い人は自分が勉強するための教材を生成AIに作ってもらって，それで練習して，レベル上げするんだと思います。例えば，
Pythonの学習用の資料を作ります。まず基礎的な文法の説明があり、その次にプログラム例を示します。プログラム例の一部を穴埋めにして、プログラムの内容を学べるようにします。
一連の内容をwebブラウザで表示して、穴埋め問題はブラウザ上でプログラムを実行して、正誤を確認できるようにします。
ってのを生成AIに投げると，ブラウザ上でPythonの勉強できる素材をすぐに作ってくれます。基礎的な勉強ができたら，応用編でプログラムつくるってのもやってみればいいと思います。お題はAIが作ってくれます。答え合わせつきで。
ちなみに学生なら，Githubで学生だという証明の手続きをするとGithubCopilotが制限なしでつかえます。