Simulating spatial PDEs in Modelica - Accessing variable values at specific times

Question

This question is somewhat related to a previous question of mine, where I didn't quite get the right solution. Link: Earlier SO-thread

I am solving PDEs which are time variant with one spatial dimension (e.g. the heat equation - see link below). I'm using the numerical method of lines, i.e. discretizing the spatial derivatives yielding a system of ODEs which are readily solved in Modelica (using the Dymola tool). My problems arise when I simulate the system, or when I plot the results, to be precise. The equations themselves appear to be solved correctly, but I want to express the spatial changes in all the discretized state variables at specific points in time rather than the individual time-varying behavior of each discrete state.

The strategy leading up to my problems is illustrated in this Youtube tutorial, which by the way is not made by me. As you can see at the very end of the tutorial, the time-varying behavior of the temperature is plotted for all the discrete points in the rod, individually. What I would like is a plot showing the temperature through the rod at a specific time, that is the temperature as a function of the spatial coordinate. My strategy to achieve this, which I'm struggling with, is: Given a state vector of N entries:

Real[N] T "Temperature";

..I would use the plotArray Dymola function as shown below.

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

Intuitively, this would yield a plot showing the temperature as a function of the spatial coordiate, or the number in the line of discrete units, to be precise. Although this command yields a result, all T-values appear to be 0 in the plot, which is definitely not the case. My question is: How can I successfully obtain and plot the temperatures at all the discrete points at a given time? Thanks in advance for your help.

The code for the problem is as indicated below.

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;

Additional edit: The simulations show clearly that for example T[3], that is the temperature of discrete segment no. 3, starts out from 30 and ends up at 70 degrees. When I write T[3] in my command window, however, I get T3 = 0.0 in return. Why is that? This is at the heart of the problem, because the plotArray function would be working if I managed to extract the actual variable values at specific times and not just 0.0.

Suggested solution: This is a rather tedious solution to achieve what I want, and I hope someone knows a better solution. When I run the simulation in Dymola, the software generates a .mat-file containing the values of the variables throughout the time of the simulation. I am able to load this file into MATLAB and manually extract the variables of my choice for plotting. For the problem above, I wrote the following command:

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

This command will plot the temperatures (as the temperatures are stored together with their derivates in the data_2 array in the .mat-file) against the respetive number of the discrete segment/element. I was really hoping to do this inside Dymola, that is avoid using MATLAB for this. For this specific problem, the amount of variables was low on account of the simplicity of this problem, but I can easily image a .mat-file which is signifanctly harder to navigate through manually like I just did.

Answer 1

Although you do not mention it explicitly I assume that you enter your plotArray command in Dymola's command window. That won't work directly, since the variables you see there do not include your simulation results: If I simulate your model, and then enter T[:] in Dymola's command window, then the printed result is

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

I'm not a Dymola expert, and the only solution I've found (to actively store and load the desired simulation results) is quite cumbersome:

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))