Gemini
Gemini
model CompressibleFDM_1D
parameter Integer N = 100 "格子点数";
parameter Real L = 1.0 "管の長さ";
parameter Real dx = L/(N-1);
parameter Real gamma = 1.4;
// --- 格子点上の保存変数 ---
Real rho[N](start={if i < N/2 then 1.0 else 0.125 for i in 1:N});
Real rho_u[N](start=fill(0.0, N));
Real E[N](start={if i < N/2 then 2.5 else 0.25 for i in 1:N});
// --- 代数変数(各点での物理量) ---
output Real p[N] "圧力";
Real u[N] "流速";
Real F[N, 3] "フラックスベクトル";
equation
// 1. 各格子点での補助変数とフラックスの計算
for i in 1:N loop
u[i] = rho_u[i] / rho[i];
p[i] = (gamma - 1) * (E[i] - 0.5 * rho[i] * u[i]^2);
// フラックスベクトルの各成分
F[i, 1] = rho[i] * u[i];
F[i, 2] = rho[i] * u[i]^2 + p[i];
F[i, 3] = (E[i] + p[i]) * u[i];
end for;
// 2. 内部格子点 (2nd to N-1) の支配方程式
for i in 2:N-1 loop
let
// 局所的な最大信号速度(数値粘性係数)
c = sqrt(gamma * p[i] / rho[i]);
alpha = abs(u[i]) + c;
in
// 質量保存
der(rho[i]) = -(F[i+1, 1] - F[i-1, 1])/(2*dx)
+ alpha * (rho[i+1] - 2*rho[i] + rho[i-1]) / (2*dx);
// 運動量保存
der(rho_u[i]) = -(F[i+1, 2] - F[i-1, 2])/(2*dx)
+ alpha * (rho_u[i+1] - 2*rho_u[i] + rho_u[i-1]) / (2*dx);
// エネルギー保存
der(E[i]) = -(F[i+1, 3] - F[i-1, 3])/(2*dx)
+ alpha * (E[i+1] - 2*E[i] + E[i-1]) / (2*dx);
end for;
// 3. 境界条件 (単純なNeumann条件: 勾配ゼロ)
// 左端 (i=1)
rho[1] = rho[2];
rho_u[1] = -rho_u[2]; // 壁条件なら速度反転、透過なら rho_u[1]=rho_u[2]
E[1] = E[2];
// 右端 (i=N)
rho[N] = rho[N-1];
rho_u[N] = rho_u[N-1];
E[N] = E[N-1];
end CompressibleFDM_1D;