Modelicaモデルを作らせてみる

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;