在 Modelica 中模拟空间偏微分方程 - 在特定时间访问变量值

Question

这个问题与我之前的问题有些相关，我没有完全得到正确的解决方案。链接：Earlier SO-thread

我正在求解具有一个空间维度的时变偏微分方程（例如热方程 - 请参见下面的链接）。我正在使用线的数值方法，即将空间导数离散化，从而产生一个在 Modelica 中很容易解决的 ODE 系统（使用 Dymola 工具）。准确地说，当我模拟系统或绘制结果时，我的问题就会出现。方程本身似乎得到了正确求解，但我想表达所有离散状态变量在特定时间点的空间变化，而不是每个离散状态的单个时变行为。

Youtube tutorial 说明了导致我的问题的策略，顺便说一下，这不是我制作的。正如您在本教程最后看到的那样，温度随时间变化的行为分别绘制在棒中的所有离散点上。我想要的是显示在特定时间通过棒的温度的图，即作为空间坐标函数的温度。我正在努力实现这一目标的策略是：给定一个包含 N 个条目的状态向量：

Real[N] T "Temperature";

..我会使用plotArray Dymola 函数，如下所示。

plotArray( {i for i in 1:N}, {T[i] for i in 1:N} )

直观地说，这将产生一个图，显示温度作为空间坐标的函数，或者准确地说是离散单位行中的数字。尽管此命令会产生结果，但图中的所有 T 值似乎都为 0，这绝对不是这种情况。我的问题是：如何在给定时间成功获取并绘制所有离散点的温度？提前感谢您的帮助。

问题的代码如下所示。

model conduction

   parameter Real rho = 1;
   parameter Real Cp = 1;
   parameter Real L = 1;
   parameter Real k = 1;
   parameter Real Tlo = 0;
   parameter Real Thi = 100;
   parameter Real Tinit = 30;
   parameter Integer N = 10 "Number of discrete segments";
   Real T[N-1] "Temperatures";
   Real deltaX = L/N;

initial equation 
   for i in 1:N-1 loop
     T[i] = Tinit;
   end for;

equation 

   rho*Cp*der(T[1]) = k*( T[2] - 2*T[1] + Thi) /deltaX^2;
   rho*Cp*der(T[N-1]) = k*( Tlo - 2*T[N-1] + T[N-2]) /deltaX^2;

   for i in 2:N-2 loop
     rho*Cp*der(T[i]) = k*( T[i+1] - 2*T[i] + T[i-1]) /deltaX^2;
   end for
   annotation (uses(Modelica(version="3.2")));

end conduction;

附加编辑： 模拟清楚地表明，例如 T[3]，即离散段的温度。 3，从30度开始，到70度结束。但是，当我在命令窗口中写入 T[3] 时，我得到 T3 = 0.0 作为回报。这是为什么？这是问题的核心，因为如果我设法在特定时间提取实际变量值而不仅仅是 0.0，plotArray 函数就会起作用。

建议的解决方案：这是实现我想要的一个相当繁琐的解决方案，我希望有人知道更好的解决方案。当我在 Dymola 中运行模拟时，软件会生成一个 .mat 文件，其中包含整个模拟期间的变量值。我能够将此文件加载到 MATLAB 中并手动提取我选择的绘图变量。针对上面的问题，我写了如下命令：

plot( [1:9]' , data_2(2:2:18 , 10)' )

此命令将绘制温度（因为温度与它们的导数一起存储在 .mat 文件中的 data_2 数组中）与离散段/元素的相应编号。我真的希望在 Dymola 中做到这一点，即避免为此使用 MATLAB。对于这个特定的问题，由于这个问题的简单性，变量的数量很少，但我可以轻松地创建一个 .mat 文件，这比我刚刚做的手动导航要困难得多。

Answer 1

虽然您没有明确提及，但我假设您在 Dymola 的命令窗口中输入了 plotArray 命令。这不会直接起作用，因为你在那里看到的变量不包括你的模拟结果：如果我模拟你的模型，然后在 Dymola 的命令窗口中输入T[:]，那么打印的结果是

T[:]
 = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0}

我不是 Dymola 专家，我发现的唯一解决方案（主动存储和加载所需的模拟结果）非常麻烦：

simulateModel("conduction", resultFile="conduction.mat")
n = readTrajectorySize("conduction.mat")
X = readTrajectory("conduction.mat", {"Time"}, n)
Y = readTrajectory("conduction.mat", {"T[1]", "T[2]", "T[3]"}, n)
plotArrays(X[1, :], transpose(Y))