state space matrix magnitude changes sampling rate

Here's the code, commented through. You don't need to know the physics to help me with this problem (although if you see a mistake I made, please let me know!). Basically, there are two problems: specifying the sample time in the ss call gives me answers I don't understand, and the matrix magnitude affects the sampling rate.

close all
clear
format long
Aoverdx = (0.005145/(2.5*10^-5));
k = 429; %thermal conductivity of silver 
kg = 0.1349; % 0.1349 kg silver
specificheatcapacity = 235; %235 J/kg K
mc = kg*specificheatcapacity; 
% The equation for current temperature difference is:
% T[t+1] =T_initial + 0*T[t] + Q/mc, where T_initial is x0(1)
% This goes in the first line of the transition matrix A.

% The equation for total heat flow is:
% Q[t+1] = 1*Q[t] + (-kA/dx)*T
% This goes in the second line of the transition matrix A.
A = [0 1/mc;-k*Aoverdx 0];
% A= A/1000;
% ^test what happens when you uncomment this line and line 26
B=[0;0];D=[0]; %not testing for inputs right now 
C = [1 0]; %observing T, the temperature difference. 
% Try switching to C = [0 1] to see how the change in heat flux 
% synchronizes with the temperature changes, displaying undamped
% energy exchange
sys1 = ss(A,B,C,D); %setting sample times gives really weird errors
t = 0:0.001:.12;
% t=t*1000;
u = zeros(size(t)); %no inputs
startingTempDiff = 20; %T_initial ,temp difference = 20 K 
x0=[startingTempDiff 0]; 
[y,t]=lsim(sys1,u,t,x0); 
lsim(sys1,u,t,x0) 
TimeToMin = t(find(y==min(y))) % notice how this time changes when 
% you uncomment lines 17 and 26

% Also notice how if you JUST uncomment line 26, an undersampling
% error appears, but if you reduce the matrix magnitude in line 17,
% that error goes away. Why does magnitude affect sampling rate? 

Ultimately, I really need the simulation to show the same TimeToMin no matter what the time vector spacing is; this is a physical process that should take a certain amount of time no matter how quickly you sample it.

Thanks for the help!

Context, if you're interested:

I'm trying to develop a controller for heat flow across a silver medium with a uniform thickness. To do that, I plan to simulate with a state space model that incorporates the physical laws, convert to transfer function, specify rise time and steady state error, and then use root locus optimization to determine the parameters I need for my controller.

Answer 1

[Jon]

I think that you are just getting a little confused about the time scaling.

In the nominal case you can see that the system reaches a minimum at t = 0.060 seconds.

In the second case (uncommenting lines 17 and 26) you divide A by 1000, since now dx/dt = A/1000*x +Bu, this makes all of your derivatives dx/dt much smaller and so this system evolves more slowly. You then counteract that by looking (sampling) at 1/1000 the sampling rate, and find that the system reaches a minimum at 60 "seconds". Your mistake is considering this second result as being in seconds, it is infact in milliseconds. (since you divided the A by 1000, and mulitplied the t by 1000).

So you actually have the same result in the second case. It reaches a minimum at t = 60 ms which is equal to 0.060 s.

For the third case, just multiplying the time by 1000, you will definitely have a problem as the time scales of the problem are such that it changes a lot in milliseconds, and you are now sampling it at intervals of 1 second.

So in general lsim doesn't know anything about what time units you are using, they could be ms, s, months, centuries, but the system matrices, A,B which define the rate of change of the state, and the time vector have to use consistent set of time units.

Usually it is best to just consider the time as being in seconds and scale your A and B matrices accordingly and then generate your time vector in seconds with a small enough increment to resolve the fastest dynamics in your plant.

[Ethan Leonard]

I appreciate the thorough answer, especially the reason why just multiplying the time scale causes relative undersampling; lsim recognizes that the speedy signal change won't be picked up by the slow sampling rate. That fully answers half my question, however:

How do you know that the minimum is at t=0.060 "seconds", specifically? I didn't specify a time interval, and my intuition tells me that if a piece of silver has a 20 degree Kelvin heat differential across a quarter centimeter thickness, it will NOT fully resolve that difference in 60 ms (actually, exactly half that to get to 0 degrees across both sides, which is the real case).

Here's a physical comparison to that situation, to explain why it doesn't seem likely. Imagine taking two thin strips of silver, one at freezing temperature and one at room temperature, and touching them together. According to this result, it would take just 30ms for them to reach the same temperature. I doubt that the time scaling is correct here.

Although it is true that the equations in the transition matrix default to seconds, ,

it doesn't seem like the result is being computed in seconds. When I tried specifying the time vector to be sampled in seconds, the result was incomprehensible (this is the second half of my question). Recomputing everything here for convenience, but all I changed was (lines 24 and 25).

close all
clear
format long
Aoverdx = (0.005145/(2.5*10^-5));
k = 429; %thermal conductivity of silver 
kg = 0.1349; % 0.1349 kg silver
specificheatcapacity = 235; %235 J/kg K
mc = kg*specificheatcapacity; 
% The equation for current temperature difference is:
% T[t+1] =T_initial + 0*T[t] + Q/mc, where T_initial is x0(1)
% This goes in the first line of the transition matrix A.

% The equation for total heat flow is:
% Q[t+1] = 1*Q[t] + (-kA/dx)*T
% This goes in the second line of the transition matrix A.
A = [0 1/mc;-k*Aoverdx 0];
% A= A/1000;
% ^test what happens when you uncomment this line and line 26
B=[1;0];D=[0]; %not testing for inputs right now 
C = [1 0]; %observing T, the temperature difference. 
% Try switching to C = [0 1] to see how the change in heat flux 
% synchronizes with the temperature changes, displaying undamped
% energy exchange
sys1 = ss(A,B,C,D,1); %setting sample times gives really weird errors
% t = 0:0.001:.12;
t = 0:1:5;
% t=t*1000;
u = zeros(size(t)); %no inputs
startingTempDiff = 20; %T_initial ,temp difference = 20 K 
x0=[startingTempDiff 0]; 
[y,t]=lsim(sys1,u,t,x0); 
lsim(sys1,u,t,x0) 

TimeToMin = t(find(y==min(y))) 

TimeToMin = 
     2

y   

y = 6×1

1.0e+08 * 0.000000200000000 0 -0.000556996987524 0 1.551228220555446 0

This is what happens when I specify a sampling time (insane magnitude!). Likewise, whatever end time I choose, the shape is the same. Also, the simulation oscillates to zero every sample, which I also don't understand. I appreciate any help, thank you for taking the time to look through this question.

[Jon]

I'm not understanding exactly what physical situation you are trying to model, and then how you are deriving your equations.

I understand that you have a sheet of silver that is undergoing some thermal transient. But what are the specifics of the situation. Here are some example situations just to give you an idea of what I mean by the specifics of the situation:

1. Sheet initially at temperature T0, both sides exposed to ambient tempertature Ta which is different than T0
2. Sheet initially at temperature T0 subjected to known heat flux q
3. As you mention in your discussion, two sheets of metal initially at temperatures T01 and T02 suddenly brought into intimate contact.

Just a quick look at your comments in your code makes it seem as if you have already discretized the system, you are describing the state T[t+1] as a function of T[t]. If so I'm not sure how you discretized them and what time step you have used for this discretization. In any case this formulation seems surprising. I would expect that you would have a set of linear ordinary differential equations which express the rate of change of the temperature dx/dt as a function of T and the inputs.

[Ethan Leonard]

Very important to be specific, so here goes. The situation is that there is a sheet initially at temperature 20C, with a fluid experiencing an endothermic reaction on one side and a fluid with temperature 20C on the other. It is believed that the fluid undergoing endothermic reaction will be at most 0C (although almost certainly it will be colder).

Now, the heat flux through a solid at any given time is dependent on the temperature gradient, , where qx is the heat flux, k is the thermal conductivity, and dT/dx is the temperature gradient. Since actionable knowledge of the temperature gradient requires sensors interspersed throughout the medium, which I wouldn't have, this equation cannot be implemented (I think). Instead, there is a simplified version which depends only on the current temperature difference across the medium.

That version (linked above), is the heat flow equation, which models the change in energy over a certain time step. Units of are in joules per second. From this, I extrapolated that the system can be discretized by , which can be easily written into a state space equation with .

Since the temperature at the infinitesimally small portion of the strip on the cold side can be said to reach 0C in zero time steps, I figured that I can model the heat flux using the conduction equation, simplified to the one dimensional heat flow equation. Likewise, the change in temperature can be modeled by the specific heat capacity of silver, , which depends on the mass m and the energy spent dQ. This is rearranged to , where mc are constants. In a way, this seems like a differential equation without the term, so that the value is fully determined by some factor, not by it's previous values. That's why this is written in the state space as .

That's how it's set up here, but please let me know if a differential equation model could be modeled well without observation, or more simply. Thanks for the thorough answer and for digging in to the physics!

[Jon]

I'm understanding your situation better now. Before getting into how to model this, I think it would help me to have even a little more context. You mentioned at the beginning that your goal is to develop a controller to control the heat flow through the silver sheet, so from the fluid at 20 deg c to the fluid experiencing an endothermic reaction. So for the control system you need to define your "process variable" the thing you are trying to hold at a setpoint, and the "manipulated variable" the thing you get to adjust.

From what you have said so far it sounds like your "process variable" is heat flux through the silver sheet. Is that correct?

So what will your manipulated variable be? That is what will you adjust to modify this heat flow?

[Ethan Leonard]

Exactly the point; this is an incomplete picture. The process variable is actually the temperature of the fluid on the other side of the silver sheet(currently at 20 C). The manipulated variable is the surface area of the sheet, which should be just large enough to absorb heat to specification. However, the simulation so far only looks at the transfer of heat across the silver sheet.

In full, the process would be:

1) Fluid undergoing endothermic reaction absorbs heat from the silver sheet (this is convection I think, , where q is the heat flux, h is the convection heat transfer coefficient, Tw is the temperature at the wall of the sheet, and Tf is the temperature of the fluid). This is not yet part of the simulation.

Heat direction: Solid boundary→Fluid

It's possible I should have simulated this step first, but since I have some reason to believe that the reaction can easily bring the silver sheet to 0C, I chose to start from the second step, considering the cold fluid to be an infinite heat sink. I'm also a little less sharp on the thermochemistry, and wanted to make sure I got the control system structure right first.

2) Heat conducts across the silver sheet.

Heat direction: Solid boundary →Solid boundary

The goal is to bring the time taken to 0 temperature difference to within specification. This is what I was aiming at in the code I wrote. So for now, I guess we can say the temperature difference across the sheet () is the process variable.

3) Silver sheet absorbs heat from the 20C fluid.

Heat direction: Fluid→Solid boundary

This is not yet part of the simulation.

Since increasing surface area A increases heat flux proportionally in the second step (using the equation , in the code written above), that's a manipulated variable that should be optimizable.

Tangential note:

In the first step (convection to the silver sheet), increasing A will also increase the quantity of reactants, which will also increase the total energy absorbed. Since I'm not clear yet on whether or not that is a linear function of A, and since I'm also not as clear on heat transfer from a solid boundary to a fluid, I didn't start there.

It is also possible to add heat directly at the final stage, which would be part of the controller. Ultimately, the temperature in the 20C fluid is the process variable of the whole system.It's maximum performance could be adjusted by adjusting the medium through which the endothermic reaction acts (the surface area of the silver). The controller can be adjusted in real time by adding reactants, and adding heat directly.

I hope that wasn't too much detail! I appreciate your hanging on to help me out with this.

[Jon]

I think you most likely just need to do a steady state (no time transient) analysis of your situation. Silver has a very high thermal conductivity, and you have only a thin sheet of silver. This means that it is likely that the temperature gradient within the silver sheet is negligible. The system will be dominated by the relatively low conduction rates through the two convective heat transfer coefficients on the surfaces of the sheet, and the relatively long time scales of the masses of fluid on either side of the sheet (unless these volumes are very small compared to the silver sheet.

You can quantify this further by looking at the Biot number for the silver sheet. This is the ratio of the thermal resistance of the sheet to the thermal resistance of the convective film. Bi = L*h/k where L is the thickness of the sheet, h is the convective heat transfer coefficient and k is the thermal conductivity of the silver sheet. When this is small, e.g. <0.1 the gradient through the sheet is neglible. The sheet temperature will just assume some intermediate value between the 20deg c fluid and the endothermic fluid, which is dictated by the h values on each side of the sheet.

Further you can ignore transient effects (in the silver sheet) if the time scales are long compared to the characteristic time for heat to diffuse through the sheet. The characteristic diffusion time for a thickness L and thermal diffusivity alpha is given by tc = L^2/alpha. So unless you're time scales of your system are very short compared to this the sheet will already have assumed a uniform throughout the thickness steady state gradient.