How to solve set of equations Matlab

Question

I have a set of equations defining a system. I have no idea how to do this in Matlab, other than the fact that it is possible, and iterations are required.

I have three equations with three unknowns (Q, T_a_out and T_b_out).

The system is divided into several sub-parts. Information about the temperature in the first and the last sub-part must be used to find the temperatures in the sub-parts in the middle. I need the Q, T_a_out and T_b_out for each sub-part.

The equations used to describe the system are as follows:

Q=U*((T_a_out-T_b_in)+(T_a_in-T_b_out))/2
Q=m_a*(cp_a_in*T_a_in-cp_a_out*T_a_out)
Q=m_b*cp_b*(T_b_out-T_b_in)

These parameters are known:

Initial T_a_in (110)
Initial T_b_in (5)
U
m_a
m_b
n (number of sub-parts)
cp_b

The cp value is a function of the temperature it belongs to:

cp_a_in is a function of the temperature: cp_a_in=function_a(T_a_in)
cp_a_out is a function of the temperature: cp_a_out=function_a(T_a_out)

The initial a value (110) has a higher value than the final b_out temperature in the last sub-part. The initial b_in value (5) has a lower value than the final a_out value.

How can I calculate the out temperatures for each sub-part in Matlab?

Answer 1

Based on your comments I've put together a quick script that should get you started in the right direction - in addition to a few notes.

The first is that what you are trying to solve is non-trivial, especially considering that the Cp(T) is highly nonlinear and on top of that diverges near the critical point.

Also, recognizing that H = CpdT, and you're using REFPROP, a more accurate calculation would be H(T_out) - H(T_in) which would incorporate temperature effects.

With that being said, by writing out your mass balances you can formulate a system of non-linear equations to solve as shown below. Since I couldn't resolve your use of UA I just substituted in an "effectiveness factor". Feel free to use your convention of choice (say, NTU) to fill that in.

I caution you that what I've posted below is incomplete. It's a starting point and there are a few issues to consider:

I am specifying n-2 temperature variables while providing n equations meaning that the bounds are likely to be satisfied as Ta_end = Tb_end (you will see what I mean if you run it and plot). This also potentially cause a discontinuity in temperature as a function of stage.

Depending on your pressures and fluids, if you are going through a 1st order phase transition you may experience some weird numerical issues.

If you can assume a constant heat capacity you can reformulate the matrix below into a much easier and solvable system. I provide this more complex example as a starting point.

function Tf = gen_mat
n = 10; % Number of stages 
A = zeros(n,2*n+2); % Coefficients matrix.

ma = 10; % Mass flow of A.
mb = 5;  % Mass flow of B.
eff = 1.0; % Effectiveness factor.
Te(1) = 383.15; % Inlet T of A.
Te(2) = 278.15; % Inlet T of B.
P = 15000; % P (kPa)

% Generate guesses 
Ti = linspace(Te(1),Te(2),n+1)';
Tg = zeros(2*n+1,1);
for i=1:n+1
    Tg(2*i-1) = Ti(i);
    Tg(2*i) = Ti(i);
end
Tg([1,2*n+1]) = []; % We are only interested in the middle sections.

% construct coefficient matrix.
for i=1:n
    A(i,2*i-1) = -ma;
    A(i,2*i) = eff*mb;
    A(i,2*i+1) = ma;
    A(i,2*i+2) = -eff*mb;
end

% solve system of nonlinear equations
Tf = fsolve(@(T)obj_fxn(T,A,P,Te), Tg, optimset('Display', 'iter'));
Tf = [Te(1);Tf;Te(2)]; % Append temperatures.
end

function b = obj_fxn(T,A,P,Te)
T = [Te(1); T; Te(2)];
for i=1:2:length(T)-1
    x(i) = refpropm('H','T',T(i),'P',P,'CO2');
    x(i+1) = refpropm('H','T',T(i+1),'P',P,'water');
end
b = A*x';
end

You can create a crude plot of the above after running it like this: plot(0:10,Tf(1:2:end)-273,10:-1:0,Tf(2:2:end)-273).

Sorry that this is not more organized, but I hope it helps.