M2の岩﨑さんの修論が日本原子力学会計算科学技術部会 部会学生論文賞に選ばれました。おめでとう。
ご挨拶
福井大学 工学部機械・システム工学科原子力安全工学コース，福井大学大学院 工学研究科安全社会基盤工学専攻原子力安全工学コース 熱水力研究室のホームページです。
当研究室では，原子炉内の熱流動に関する解析的・実験的研究を行っています。
- 解析ではBest Estimateコードと呼ばれる，原子炉を1次元の流路で模擬したコード及び解析対象を詳細なメッシュで分割して解析を行う数値流体力学（CFD）コードを用いて研究を進めています。
- データ同化と呼ばれる観測データ（実験データ）と解析データの双方を用いて，より精度の高い解析結果を求める手法の研究も行っています。
- 高温ガス炉や核融合炉関係の熱流動の研究も進めようとしています。
詳細については右上のRESEARCHタブもしくはPAPERSタブをクリックしてください。
学生募集
修士課程・博士課程の大学院生を募集しています。原子炉の熱流動に興味がある方はぜひご連絡ください。
研究室見学
研究室見学はいつでも歓迎します。 連絡は ishigaki@u-fukui (後ろに.ac.jpをつけて)までよろしくお願いします。
生成AIによるプログラム生成
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]春の原子力学会
日本原子力学会2026年春の年会
熊本の原子力学会に来ました。去年の秋に続けて，九州です。 M2の岩﨑さんが修了間際にかかわらず，データ同化による乱流シュミット数の推定について発表してくれました。
石垣は企画セッションの座長のお仕事をしました。AI関連の企画セッションでしたが，AI絡みの企画セッションはいくつか出てました。生成AIがこれだけ使えるようになると，当然かという気もします。
熊本駅前はきれいでした。ご飯も美味しかったのでまた来たいです。いたるところにくまモン。
写真にはないですが，いつにもない数の企業ブースの数でした。
NSRRのプロジェクションマッピング用のモデルを作ったそうです。JAEA時代の同期が紹介してました。
2025年度卒論発表・修論発表
卒論発表・修論発表
今週，卒論発表と修論発表がありました。
無事終わりました。みなさん，お疲れ様でした。
修論発表で岩﨑さんが優秀論文賞を受賞しました。おめでとうございます。
表彰式の様子
CFD4NRS-10 @水戸
国際会議CFD4NRS−10
原子力のCFDの会議CFD4NRS-10 (CFD for Nuclear Reactor Safety)が水戸でありました。はじめての日本開催。
運営はJAEAの安全研究センター熱水力安全研究グループのみなさま。準備と運営お疲れ様でした。
凝縮実験のデータ同化解析について石垣が発表しました。
会場の水戸市民会館からは水戸のクネクネタワーが近くに見えました。
ベトナムでの集中講義
ハノイ工科大学での集中講義
ツイニングプログラムというのがあり，ベトナムのハノイ工科大学から日本の参加大学に編入できるというシステムです。これに福井大も参加しており，集中講義に行ってきました。
担当は「熱力学」でした。日本に留学するのが目的のプログラムなので，学生さんは日本語の勉強をしていて，講義は日本語です。 月水金は4コマ（1コマ45分），火木は3コマで，合計5日間のわりと長めの集中講義です。
今週ずっとベトナムですが，先週の木金は筑波（製作中の実験装置の確認）にいってからのベトナムなので，結構大変です。。。
再来週はCFD4NRSで水戸です。。。
相変わらずのカオスな道路事情
大学構内
ベトナムのKFC。パスタとかまで売ってます。オリジナルチキンはいまいちでした。
2025年の原子力学会
日本原子力学会2025年秋の大会
小倉の原子力学会に来ました。 M1の和田さんがデータ同化の発表しました。最終日の最終セッションのわりに人がそれなりにいてくれて良かったです。
看板に予算を使ったみたいです。
北流研
第75回北陸流体工学研究会
金沢工業大学で第75回北陸流体工学研究会があって，データ同化について発表してきました。 北陸流体工学研究会は2回目なので，自己紹介的なつもりの発表でした。
初めての金沢。北陸新幹線速すぎて，敦賀から金沢で1時間かからないのは衝撃でした。
Volume Penalization法のテスト2
テストその2 角柱後流のカルマン渦
今回は角柱後流のカルマン渦の計算をVP法と通常の計算とで実施して，渦の放出周期が合うかやってみました。 完全に四角なので，一致するはず。
レイノルズ数は100です。
通常の計算（角柱部分のメッシュを設定）による流速分布
VP法による計算による流速分布
VP法のほうが不安定化（流れが振動しはじめて，渦が出るまで）が早いですが，最終的な渦の放出周期はほぼ一致しています。 （ここまで来るのに多少手間取ったりもしましたが。。。）
下図はスパン方向速度の時系列データです。赤が通常の計算，黒がVP法の計算によるものです。
ついでに渦度の動画
通常の計算（角柱部分のメッシュを設定）
VP法による計算
Volume Penalization法のテスト
VP法のテスト
境界埋め込み法の1つであるVP法のコードをOpenFOAMに実装しました。実装したコードはこちら。 https://qiita.com/ishigaki/items/cfb20bc58c2ae6a633a0
境界適合格子でやったときとちゃんと比べてみました。NACA0012翼周りの流れをやります。
条件
- 流入速度44m/s
- 迎角4度（境界速度の角度を調整して迎角を設定）
- コード長は1m
- Reynolds数は粘性係数で調整。Re=UinC/nu。Uは流入速度，Cはコード長，nuは動粘性係数。
メッシュ (blockMeshDict)
pisoFoam用のメッシュと，VP法つかったpisoFaomのvppisoFoam用のメッシュを作りました。blockMeshだけでほぼNACA0012通りの形状を再現できます。
メッシュ幅はおおよそ5mmです。
境界適合格子（ここをクリックするとblockMeshDictの詳細が出ます）
polyLineで翼表面の座標を直線で補間します。
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 2.3.1 |
| \\ / A nd | Web: www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
FoamFile
{
version 2.0;
format ascii;
class dictionary;
object blockMeshDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
convertToMeters 1;
vertices
(
( -0.5 -1.25 -1 )
( -0.5 1.25 -1 )
( 0.3 -0.06001727 -1 )
( 0.3 0.06001727 -1 )
( 1.0089 -1.25 -1 )
( 1.0089 0 -1 )
( 1.0089 1.25 -1 )
( 3 -1.25 -1 )
( 3 0 -1 )
( 3 1.25 -1 )
( -0.5 -1.25 1 )
( -0.5 1.25 1 )
( 0.3 -0.06001727 1 )
( 0.3 0.06001727 1 )
( 1.0089 -1.25 1 )
( 1.0089 0 1 )
( 1.0089 1.25 1 )
( 3 -1.25 1 )
( 3 0 1 )
( 3 1.25 1 )
);
blocks
(
hex (0 2 3 1 10 12 13 11) (250 120 1) simpleGrading (1 1 1)
hex (0 4 5 2 10 14 15 12) (142 250 1) simpleGrading (1 1 1)
hex (3 5 6 1 13 15 16 11) (142 250 1) simpleGrading (1 1 1)
hex (4 7 8 5 14 17 18 15) (398 250 1) simpleGrading (1 1 1)
hex (5 8 9 6 15 18 19 16) (398 250 1) simpleGrading (1 1 1)
);
edges
(
polyLine 2 3 (
( 0.28 -0.05992581 -1 )
( 0.26 -0.059636853 -1 )
( 0.24 -0.059131171 -1 )
( 0.22 -0.058386276 -1 )
( 0.2 -0.05737543 -1 )
( 0.18 -0.056066189 -1 )
( 0.16 -0.054418208 -1 )
( 0.14 -0.052379745 -1 )
( 0.12 -0.049881816 -1 )
( 0.1 -0.046827704 -1 )
( 0.08 -0.043072299 -1 )
( 0.06 -0.03837581 -1 )
( 0.04 -0.032277225 -1 )
( 0.02 -0.023597771 -1 )
( 0.01 -0.017037074 -1 )
( 0.005 -0.012213147 -1 )
( 0 0 -1 )
( 0.005 0.012213147 -1 )
( 0.01 0.017037074 -1 )
( 0.02 0.023597771 -1 )
( 0.04 0.032277225 -1 )
( 0.06 0.03837581 -1 )
( 0.08 0.043072299 -1 )
( 0.1 0.046827704 -1 )
( 0.12 0.049881816 -1 )
( 0.14 0.052379745 -1 )
( 0.16 0.054418208 -1 )
( 0.18 0.056066189 -1 )
( 0.2 0.05737543 -1 )
( 0.22 0.058386276 -1 )
( 0.24 0.059131171 -1 )
( 0.26 0.059636853 -1 )
( 0.28 0.05992581 -1 )
)
polyLine 12 13 (
( 0.28 -0.05992581 1 )
( 0.26 -0.059636853 1 )
( 0.24 -0.059131171 1 )
( 0.22 -0.058386276 1 )
( 0.2 -0.05737543 1 )
( 0.18 -0.056066189 1 )
( 0.16 -0.054418208 1 )
( 0.14 -0.052379745 1 )
( 0.12 -0.049881816 1 )
( 0.1 -0.046827704 1 )
( 0.08 -0.043072299 1 )
( 0.06 -0.03837581 1 )
( 0.04 -0.032277225 1 )
( 0.02 -0.023597771 1 )
( 0.01 -0.017037074 1 )
( 0.005 -0.012213147 1 )
( 0 0 1 )
( 0.005 0.012213147 1 )
( 0.01 0.017037074 1 )
( 0.02 0.023597771 1 )
( 0.04 0.032277225 1 )
( 0.06 0.03837581 1 )
( 0.08 0.043072299 1 )
( 0.1 0.046827704 1 )
( 0.12 0.049881816 1 )
( 0.14 0.052379745 1 )
( 0.16 0.054418208 1 )
( 0.18 0.056066189 1 )
( 0.2 0.05737543 1 )
( 0.22 0.058386276 1 )
( 0.24 0.059131171 1 )
( 0.26 0.059636853 1 )
( 0.28 0.05992581 1 )
)
polyLine 3 5 (
( 0.32 0.05992788 -1 )
( 0.34 0.059672249 -1 )
( 0.36 0.059263278 -1 )
( 0.38 0.058712465 -1 )
( 0.4 0.058030108 -1 )
( 0.42 0.057225479 -1 )
( 0.44 0.056306948 -1 )
( 0.46 0.055282094 -1 )
( 0.48 0.054157786 -1 )
( 0.5 0.052940252 -1 )
( 0.52 0.051635135 -1 )
( 0.54 0.050247543 -1 )
( 0.56 0.048782083 -1 )
( 0.58 0.047242897 -1 )
( 0.6 0.045633691 -1 )
( 0.62 0.043957754 -1 )
( 0.64 0.042217983 -1 )
( 0.66 0.040416898 -1 )
( 0.68 0.038556656 -1 )
( 0.7 0.036639067 -1 )
( 0.72 0.034665604 -1 )
( 0.74 0.032637411 -1 )
( 0.76 0.030555316 -1 )
( 0.78 0.028419835 -1 )
( 0.8 0.02623118 -1 )
( 0.82 0.023989265 -1 )
( 0.84 0.021693713 -1 )
( 0.86 0.019343859 -1 )
( 0.88 0.016938755 -1 )
( 0.9 0.014477173 -1 )
( 0.92 0.011957609 -1 )
( 0.94 0.009378289 -1 )
( 0.96 0.006737166 -1 )
( 0.98 0.004031929 -1 )
( 1 0.00126 -1 )
)
polyLine 2 5 (
( 0.32 -0.05992788 -1 )
( 0.34 -0.059672249 -1 )
( 0.36 -0.059263278 -1 )
( 0.38 -0.058712465 -1 )
( 0.4 -0.058030108 -1 )
( 0.42 -0.057225479 -1 )
( 0.44 -0.056306948 -1 )
( 0.46 -0.055282094 -1 )
( 0.48 -0.054157786 -1 )
( 0.5 -0.052940252 -1 )
( 0.52 -0.051635135 -1 )
( 0.54 -0.050247543 -1 )
( 0.56 -0.048782083 -1 )
( 0.58 -0.047242897 -1 )
( 0.6 -0.045633691 -1 )
( 0.62 -0.043957754 -1 )
( 0.64 -0.042217983 -1 )
( 0.66 -0.040416898 -1 )
( 0.68 -0.038556656 -1 )
( 0.7 -0.036639067 -1 )
( 0.72 -0.034665604 -1 )
( 0.74 -0.032637411 -1 )
( 0.76 -0.030555316 -1 )
( 0.78 -0.028419835 -1 )
( 0.8 -0.02623118 -1 )
( 0.82 -0.023989265 -1 )
( 0.84 -0.021693713 -1 )
( 0.86 -0.019343859 -1 )
( 0.88 -0.016938755 -1 )
( 0.9 -0.014477173 -1 )
( 0.92 -0.011957609 -1 )
( 0.94 -0.009378289 -1 )
( 0.96 -0.006737166 -1 )
( 0.98 -0.004031929 -1 )
( 1 -0.00126 -1 )
)
polyLine 13 15 (
( 0.32 0.05992788 1 )
( 0.34 0.059672249 1 )
( 0.36 0.059263278 1 )
( 0.38 0.058712465 1 )
( 0.4 0.058030108 1 )
( 0.42 0.057225479 1 )
( 0.44 0.056306948 1 )
( 0.46 0.055282094 1 )
( 0.48 0.054157786 1 )
( 0.5 0.052940252 1 )
( 0.52 0.051635135 1 )
( 0.54 0.050247543 1 )
( 0.56 0.048782083 1 )
( 0.58 0.047242897 1 )
( 0.6 0.045633691 1 )
( 0.62 0.043957754 1 )
( 0.64 0.042217983 1 )
( 0.66 0.040416898 1 )
( 0.68 0.038556656 1 )
( 0.7 0.036639067 1 )
( 0.72 0.034665604 1 )
( 0.74 0.032637411 1 )
( 0.76 0.030555316 1 )
( 0.78 0.028419835 1 )
( 0.8 0.02623118 1 )
( 0.82 0.023989265 1 )
( 0.84 0.021693713 1 )
( 0.86 0.019343859 1 )
( 0.88 0.016938755 1 )
( 0.9 0.014477173 1 )
( 0.92 0.011957609 1 )
( 0.94 0.009378289 1 )
( 0.96 0.006737166 1 )
( 0.98 0.004031929 1 )
( 1 0.00126 1 )
)
polyLine 12 15 (
( 0.32 -0.05992788 1 )
( 0.34 -0.059672249 1 )
( 0.36 -0.059263278 1 )
( 0.38 -0.058712465 1 )
( 0.4 -0.058030108 1 )
( 0.42 -0.057225479 1 )
( 0.44 -0.056306948 1 )
( 0.46 -0.055282094 1 )
( 0.48 -0.054157786 1 )
( 0.5 -0.052940252 1 )
( 0.52 -0.051635135 1 )
( 0.54 -0.050247543 1 )
( 0.56 -0.048782083 1 )
( 0.58 -0.047242897 1 )
( 0.6 -0.045633691 1 )
( 0.62 -0.043957754 1 )
( 0.64 -0.042217983 1 )
( 0.66 -0.040416898 1 )
( 0.68 -0.038556656 1 )
( 0.7 -0.036639067 1 )
( 0.72 -0.034665604 1 )
( 0.74 -0.032637411 1 )
( 0.76 -0.030555316 1 )
( 0.78 -0.028419835 1 )
( 0.8 -0.02623118 1 )
( 0.82 -0.023989265 1 )
( 0.84 -0.021693713 1 )
( 0.86 -0.019343859 1 )
( 0.88 -0.016938755 1 )
( 0.9 -0.014477173 1 )
( 0.92 -0.011957609 1 )
( 0.94 -0.009378289 1 )
( 0.96 -0.006737166 1 )
( 0.98 -0.004031929 1 )
( 1 -0.00126 1 )
)
);
boundary
(
inlet
{
type patch;
faces
(
(0 1 11 10)
);
}
outlet
{
type patch;
faces
(
(0 4 14 10)
( 4 7 17 14)
(1 6 16 11)
(6 9 19 16)
(7 8 18 17)
(8 9 19 18)
);
}
airfoil
{
type wall;
faces
(
(2 5 15 12)
(2 3 13 12)
(3 5 15 13)
);
}
front
{
type empty;
faces
(
(0 2 3 1)
(0 4 5 2)
(3 5 6 1)
(4 7 8 5)
(5 8 9 6)
(10 12 13 11)
(10 14 15 12)
(13 15 16 11)
(14 17 18 15)
(15 18 19 16)
);
}
);
mergePatchPairs
(
);
// ************************************************************************* //