Vary a parameter and plot the output

Hi, I have the following script for which I am trying to determine how the outlet temperature (Tout) varies depending on the solar radiation (I) available. I have tried to set (I) as an array and run the code as a for loop but this doesn't seem to give the desired result. Additionally it doesn't plot (I) against (Tout) at the end. Can anyone provide me with some insight as to how to rectify this. I have attached the code.

Kind regards

Answer 1

[Star Strider]

Try this:

function main
x0 = [500 500 500 500 5 5];
opt = optimset('Display','off');  % 'display', 'iter'
I = [500 550 600 650 700 750]; % Solar radiation
for k = 1:numel(I)
    res(k,:) = fsolve(@(x)f(x,I(k)),x0,opt);
%     res'
end

figure
plot(I, res(:,3))
grid
xlabel('Solar Radiation')
ylabel('T_{out}')
end

That ran for me without error in R2020a.

[kasim]

Hi, where would you place the for loop as it seems that you have moved it into the first function?

[Star Strider]

I forgot to include the slightly changed function definition, that being:

function res = f(x,I)

Then eliminate the definition of ‘I’ within the fuction, since it would otherwise overwrite the value of the passed parameter ‘I’. (The rest of the function code is unchanged.)

I moved the loop to the main function because it belongs there. The for loop uses every element of ‘I’, once in each iteration, passing it to the function as an additional parameter.

The full code is then:

function main
x0 = [500 500 500 500 5 5];
opt = optimset('Display','off');  % 'display', 'iter'
I = [500 550 600 650 700 750]; % Solar radiation
for k = 1:numel(I)
    res(k,:) = fsolve(@(x)f(x,I(k)),x0,opt);
%     res'
end

figure
plot(I, res(:,3))
grid
xlabel('Solar Radiation')
ylabel('T_{out}')
end

function res = f(x,I)
[Tg Tp Tout Tm hrpg hrgs] = deal(x(1),x(2),x(3),x(4),x(5),x(6));
 
%for I = [500 550 600 650 700 750]; % Solar radiation
% I = 750;
Ta = 295;        % ambient temperature
Tin = Ta;        % inlet temperature
m = 0.02;        % mass flowrate
cp = 1009;          % specific heat capacity of air
ag = 0.01;       % glass absorption
taug = 0.91;     % glass transmissivity
ap = 0.95;       % plate absorption
Eg = 0.96;       % glass emittance
Ep = 0.3;        % plate emittance
Sigma = 5.67E-8; % Boltzman constant
L = 1.5;         % Collector Length
W = 1.2;         % Collector Width
A=L*W;           % Collector Area
hw = 5;
hp = 10;
hg = 10;
 
eq1 = ag*I + hg*(Tm-Tg) + hrpg*(Tp-Tg)- ((hw+hrgs)*(Tg-Ta));
eq2 = m*cp*(Tout-Tin)-(hg*A*(Tg-Tm))-(hp*A*(Tp-Tm));
eq3 = hp*(Tp-Tm)+hrpg*(Tp-Tg)-ap*taug*I;
eq4 = 2*Tm-Tin-Tout;
eq5 = hrpg - (Sigma*((Tp^2)+(Tg^2))*(Tp+Tg))/((1/Ep)+(1/Eg)-1);
eq6 = hrgs - Sigma*Eg*((Tg^4)-(Ta^4))/(Tg-Ta) ;
res = [eq1 eq2 eq3 eq4 eq5 eq6];
% vars = symvar(eqs);
% S = solve(eqs,vars);
end

[kasim]

Thank you very much! It works just as needed.

[Star Strider]

As always, my pleasure!

[kasim]

Sorry one last question, if I now want to make further calculations (i.e. determine efficiency) using the Tout that has just been calculated. Where should I place that in the code? Again, thank you!

[Star Strider]

As always, my pleasure!

Since ‘Tout’ is a vector created in the loop, put all calculations involving ‘Tout’ anyplace after the loop, similar to where the plot call is now. (Those calculations can come before or after the plot call, depending on what you want to do with them. They must be between the loop and the function main end statement.) Another loop is likely not necessary, since once ‘Tout’ is created, most calculations involving it can be done with vectorized code. .

[kasim]

Yes that was my initial thought too. When I attempt to add a calculation (i.e efficiency) before the main function ends as below: I get the error shown:

because those variables are stored in the main body.

Essentially for each iteration of the solar radiation, I would like to compute its efficiency against the outlet temperature. So would I have to add it somewhere in the main body.

Regards

n = m*cp*(Tout-Tin)/I*A


Undefined function or variable 'm'.

[Star Strider]

Copy ‘m’, ‘cp’, ‘A’ and ‘Ta’ to the main function (since they don’t change) and then do the calculation. Only ‘Tout’ and ‘I’ are vectors, so with those changes and using element-wise division, the ‘n’ calculation should work as written.

Note that the expression divides by ‘I’ (and then multiplies the numerator by ‘A’), so element-wise division willl be necessary, both because of that and because ‘Tout’ and ‘I’ are same-sized vectors:

n = m*cp*(Tout-Tin)./I*A;

I didn’t test that, however it should work.

[kasim]

Tout = res(:,3)

%To calculate the efficiency

m = 0.02;

cp = 1009;

Tin = 295;

A = 1.8;

n = m*cp*(Tout-Tin)./I*A

For some reason it doesn't work, MATLAB calculates n for each value in the matrix and outputs that rather than just calculating it for Tout.

[Star Strider]

The ‘Tout’ vector have values for every value of ‘I’, so ‘n’ will be a vector that represents every value of ‘Tout./I’, as calculated in the loop, scaled by the constants.

That was the entire point of calculating ‘Tout’ as a function of ‘I’.

[kasim]

n =

    2.0508    1.8644    1.7090    1.5775    1.4649    1.3672
    2.2534    2.0486    1.8778    1.7334    1.6096    1.5023
    2.4556    2.2323    2.0463    1.8889    1.7540    1.6370
    2.6573    2.4157    2.2144    2.0440    1.8980    1.7715
    2.8585    2.5986    2.3821    2.1988    2.0418    1.9056
    3.0592    2.7811    2.5493    2.3532    2.1852    2.0395

  

Apologies, what I was trying to say was that it calculates n for not only Tout, but also Tm,Tp,Tg etc as can be seen by the following code. Is there a way to fix this so it only gives the values corresponding to Tout? i.e. row 3

[Star Strider]

The ‘Tout’ value is ‘res(:,3)’. Define it that way (perhaps as ‘Toutv’) and you should get the result you want:

Toutv = res(:,3);
n = m*cp*(Toutv-Tin)./I*A;

That should produce a vector rather than a matrix.