How do you solve a System of 11 Non-Linear Equations Inter-linked to one another?

Dear All,

I am currently trying to solve a system of non-linear equations to calculate the overall heat transfer coefficient of a triple concentric tube heat exchanger. Due to the way the equations are inter-linked to one another, I am unsure as to where to start my code with.

Problem Definition:

* U1 and U2 are the two variables I would like to solve for.

* The variables delta_T1(0), delta_T1(L), delta_T2(0), delta_T2(L), A_lm1, A_lm2, C_h, C_c1, C_c2, P1, P2, and L are known values.

The system of equations I would like to solve are as follows:

where f(L) and g(L) are defined as:

where G1, G2, G3, and G4 are solved by the following simultaneous equations:

where lambda_1 and lambda_2 are defined as:

and lastly, the variables B and C are defined as:

Attached below are the values of the known values I outlined, which hopefully would help describe the problem better:

% Testing for Parallel Flow, flow rate = 0.09kg/s

T_h_0 = 343;
T_h_L = 327;
T_c1_0 = 283;
T_c1_L = 290;
T_c2_0 = 291;
T_c2_L = 299;

r1_o = (13.51e-3)/2;
r1_i = (7.97e-3)/2;
r2_o = (45.26e-3)/2;
r2_i = (39.72e-3)/2;

% delta_T1 and delta_T2 calcs
delta_T1_0 = T_h_0 - T_c1_0;
delta_T1_L = T_h_L - T_c1_L;
delta_T2_0 = T_h_0 - T_c2_0;
delta_T2_L = T_h_L - T_c2_L;

% P1 and P2 calcs
P1 = 2 * pi * (r1_o - r1_i) / log(r1_o/r1_i);
P2 = 2 * pi * (r2_o - r2_i) / log(r2_o/r2_i);

% A_lm1 and A_lm2 calcs
A_lm1 = P1 * 2.5;
A_lm2 = P2 * 2.5;

% C_c1, C_c2, and C_h value calcs = mass flow rate x specific heat
C_c1 = 0.1 * 4200;
C_c2 = 0.1 * 4180;
C_h = 0.09 * 4190;

Any suggestions on how I could resolve this would be very much appreciated. Thank you in advance!

[Torsten]

What is L ?

Any reasonable initial guesses for U1 and U2 ?

[Christiantyo Moeprodjo]

Hi,

Thank you for the swift response!

L is the length of the heat exchanger which is 2.5 metres.

From the simulation I ran;

U1 = 663.246 ≈ 663

U2 = 255.375 ≈ 255

If I was to take a reasonable guess for U1 and U2, it would be them!

If you have further questions on any parameters that seem to be unclear, please let me know and I will try my best to clarify them for you. Thank you once again for the help!

Answer 1

[Alan Stevens]

Here's a possible logic (though I seem to find negative values for U):

U1 = 500;  U2 = 500;  % Initial guesses
U0 = [U1 U2];

[U, fval] = fminsearch(@eqns, U0);
U1 = U(1); U2 = U(2);


disp([U1 U2])

   1.0e+03 *

   -2.1221   -0.6021

disp(fval)

   6.3809e-09

function F = eqns(U)

        % Testing for Parallel Flow, flow rate = 0.09kg/s

        L = 2.5;

        T_h_0 = 343;
        T_h_L = 327;
        T_c1_0 = 283;
        T_c1_L = 290;
        T_c2_0 = 291;
        T_c2_L = 299;

        r1_o = (13.51e-3)/2;
        r1_i = (7.97e-3)/2;
        r2_o = (45.26e-3)/2;
        r2_i = (39.72e-3)/2;

        % delta_T1 and delta_T2 calcs
        delta_T1_0 = T_h_0 - T_c1_0;
        delta_T1_L = T_h_L - T_c1_L;
        delta_T2_0 = T_h_0 - T_c2_0;
        delta_T2_L = T_h_L - T_c2_L;

        % P1 and P2 calcs
        P1 = 2 * pi * (r1_o - r1_i) / log(r1_o/r1_i);
        P2 = 2 * pi * (r2_o - r2_i) / log(r2_o/r2_i);

        % A_lm1 and A_lm2 calcs
        A_lm1 = P1 * 2.5;
        A_lm2 = P2 * 2.5;

        % C_c1, C_c2, and C_h value calcs = mass flow rate x specific heat
        C_c1 = 0.1 * 4200;
        C_c2 = 0.1 * 4180;
        C_h = 0.09 * 4190;


         U1 = U(1);  U2 = U(2);
         
         B = U1*P1*(1/C_c1 - 1/C_h) + U2*P2*(1/C_c2 - 1/C_h);
         C = U1*P1*U2*P2*(1/(C_c1*C_c2) - 1/(C_h*C_c1) - 1/(C_h*C_c2));
         
         d = sqrt(B^2 - 4*C);
         lambda1 = (-B + d)/2;  lambda2 = (-B - d)/2;
         
         M = [1 1; exp(lambda1*L) exp(lambda2*L)];
         V1 = [delta_T1_0; delta_T1_L];
         V2 = [delta_T2_0; delta_T2_L];
         
         G12 = M\V1;  G1 = G12(1); G2 = G12(2);
         G34 = M\V2;  G3 = G34(1); G4 = G34(2);
         
         logratio1 = log((G1*exp(lambda1*L)+G2*exp(lambda2*L))/(G1+G2));
         logratio2 = log((G3*exp(lambda1*L)+G4*exp(lambda2*L))/(G3+G4));
         f = (L*(G1*G4*lambda1-G2*G3*lambda2) + (G2*G3-G1*G4)*logratio1)...
             /(G1*G2*(lambda1-lambda2));
         g = (L*(G2*G3*lambda1-G1*G4*lambda2) + (G1*G4-G2*G3)*logratio2)...
             /(G3*G4*(lambda1-lambda2));
         
         lndT1 = U1*A_lm1*(1/C_h - 1/C_c1) + U2*P2*f/C_h;
         lndT2 = U2*A_lm2*(1/C_h - 1/C_c2) + U1*P1*g/C_h;
         
         % lndT1 should be the same as logratio1 and
         % lndT2 should be the same as logratio2, so
         % calculate the least squares difference between 
         % each and sum to get F, which one would like to be zero.
         
         F = norm(lndT1-logratio1) + norm(lndT2-logratio2);
         
end

Answer 2

[Torsten]

U0 = [663 255];
options = optimset('TolX',1e-12,'TolFun',1e-12);
format long
U = fsolve(@fun,U0,options)

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

U = 
   1.0e+03 *

 -2.122106380642823 - 0.000000000000000i -0.602069349598804 - 0.000000000000000i

fun(U)

ans = 
   1.0e-12 *

  0.350719453479087 + 0.000000269903231i  0.611843908870924 + 0.000000646844817i

function res = fun(U)
L = 2.5;
% Testing for Parallel Flow, flow rate = 0.09kg/s
T_h_0 = 343;
T_h_L = 327;
T_c1_0 = 283;
T_c1_L = 290;
T_c2_0 = 291;
T_c2_L = 299;
r1_o = (13.51e-3)/2;
r1_i = (7.97e-3)/2;
r2_o = (45.26e-3)/2;
r2_i = (39.72e-3)/2;
% delta_T1 and delta_T2 calcs
delta_T1_0 = T_h_0 - T_c1_0;
delta_T1_L = T_h_L - T_c1_L;
delta_T2_0 = T_h_0 - T_c2_0;
delta_T2_L = T_h_L - T_c2_L;
% P1 and P2 calcs
P1 = 2 * pi * (r1_o - r1_i) / log(r1_o/r1_i);
P2 = 2 * pi * (r2_o - r2_i) / log(r2_o/r2_i);
% A_lm1 and A_lm2 calcs
A_lm1 = P1 * 2.5;
A_lm2 = P2 * 2.5;
% C_c1, C_c2, and C_h value calcs = mass flow rate x specific heat
C_c1 = 0.1 * 4200;
C_c2 = 0.1 * 4180;
C_h = 0.09 * 4190;
U1 = U(1);
U2 = U(2);
B = U1*P1*(1/C_c1-1/C_h) + U2*P2*(1/C_c2-1/C_h);
C = U1*P1*U2*P2*(1/(C_c1*C_c2)-1/(C_h*C_c1)-1/(C_h*C_c2));
lambda1 = (-B+sqrt(B^2-4*C))/2;
lambda2 = (-B-sqrt(B^2-4*C))/2;
Amat = [1 1 0 0;exp(lambda1*L) exp(lambda2*L) 0 0;0 0 1 1 ;0 0 exp(lambda1*L) exp(lambda2*L)];
bmat = [delta_T1_0;delta_T1_L;delta_T2_0;delta_T2_L];
G = Amat\bmat;
G1 = G(1);
G2 = G(2);
G3 = G(3);
G4 = G(4);
fL = (L*(G1*G4*lambda1 - G2*G3*lambda2) + (G2*G3 - G1*G4)*log((G1*exp(lambda1*L)+G2*exp(lambda2*L))/(G1+G2)))/(G1*G2*(lambda1-lambda2));
gL = (L*(G2*G3*lambda1 - G1*G4*lambda2) + (G1*G4 - G2*G3)*log((G3*exp(lambda1*L)+G4*exp(lambda2*L))/(G3+G4)))/(G3*G4*(lambda1-lambda2));
res(1) = log(delta_T1_L/delta_T1_0) - (U1*A_lm1*(1/C_h-1/C_c1) + U2*P2/C_h*fL);
res(2) = log(delta_T2_L/delta_T2_0) - (U2*A_lm2*(1/C_h-1/C_c2) + U1*P1/C_h*gL);
end

[Torsten]

@Christiantyo Moeprodjo comment moved here:

Hi,

Apologies for getting back to this quite late. However, I am having issues with the values the script is generating.

The initial guesses I took for U1 and U2 are based on a CFD simulation study I conducted for the pipe flow I'm interested in. Hence, the final values I should be getting from these set of equations should be similar to the initial guesses. However, the final values for U generated by the script are somewhat large (~3000 for U1, and ~900 for U2) compared to what I am expecting (~630 for U1, and ~450 for U20).

Do you have any idea as to why this might be the case?

Furthermore, the script outputs values for 'U' and 'ans'. What does 'ans' represent?

Thank you very much for the help, I really appreciate it.

Kind regards,

Chris

[Torsten]

Furthermore, the script outputs values for 'U' and 'ans'. What does 'ans' represent?

If you print out the value of an expression that has no explicit name, MATLAB uses ans(wer) as name for this expression:

36*25

ans = 900

a = 36*25

a = 900

[Christiantyo Moeprodjo]

That makes so much more sense now, thanks!

However, do you have any clue why the 'U' values the script generates are extremely huge compared to the initial guesses taken?

[Torsten]

However, do you have any clue why the 'U' values the script generates are extremely huge compared to the initial guesses taken?

Not really. My guess is that the function is very flat and large changes in U are necessary to make it go to zero.

But you have an option: Plot the function in the region of interest using a surface plot (i.e. plot fun(U(1),U(2)) against U(1) and U(2)). Then you can see how it behaves.

Look at the examples provided under

https://de.mathworks.com/help/matlab/ref/surf.html

to understand how to proceed.

[Christiantyo Moeprodjo]

Also, I missed a huge part of this inquiry. I have amended the delta_T relationships, which should replicate my scenario of interest better.

The updated code are as follows:

% Overall Heat Transfer Coefficient [U1, U2] in a Triple Tube Heat Exchanger

U0 = [663 255];
fun(U0)

ans = 1×2

0.5638 0.7867

options = optimset('TolX',1e-12,'TolFun',1e-12);
format long
U = fsolve(@fun,U0,options)

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

U = 1×2

1.0e+03 * 3.296358855124808 0.971517845099225

fun(U)

ans = 1×2

1.0e-11 * 0.448985293388660 -0.950273193467410

function res = fun(U)
L = 2.5;

% Testing for Counter Flow, Re = 1634
% Temperature Values [degC]
T_h_0 = 70;
T_h_L = 35.71505;
T_c1_0 = 10;
T_c1_L = 15.0687;
T_c2_0 = 18;
T_c2_L = 22.29458;

% radius of each wall (inner and outer of respective tubes)
r1_o = (13.51e-3)/2;
r1_i = (7.97e-3)/2;
r2_o = (45.26e-3)/2;
r2_i = (39.72e-3)/2;

% delta_T1 and delta_T2 calcs
delta_T1_0 = T_h_L - T_c1_0; % outlet of hot fluid - inlet of cold fluid
delta_T1_L = T_h_0 - T_c1_L; % inlet of hot fluid - outlet of cold fluid
delta_T2_0 = T_h_L - T_c2_0; % outlet of hot fluid - inlet of normal fluid
delta_T2_L = T_h_0 - T_c2_L; % inlet of hot fluid - outlet of normal fluid

% P1 and P2 calcs
P1 = 2 * pi * (r1_o - r1_i) / log(r1_o/r1_i);
P2 = 2 * pi * (r2_o - r2_i) / log(r2_o/r2_i);

% A_lm1 and A_lm2 calcs
A_lm1 = P1 * 2.5;
A_lm2 = P2 * 2.5;

% C_c1, C_c2, and C_h value calcs = mass flow rate x specific heat
C_c1 = 0.1 * 4200;
C_c2 = 0.1 * 4180;
C_h = 0.09 * 4190;

U1 = U(1);
U2 = U(2);

B = U1*P1*(1/C_c1-1/C_h) + U2*P2*(1/C_c2-1/C_h);
C = U1*P1*U2*P2*(1/(C_c1*C_c2)-1/(C_h*C_c1)-1/(C_h*C_c2));

lambda1 = (-B+sqrt(B^2-4*C))/2;
lambda2 = (-B-sqrt(B^2-4*C))/2;

Amat = [1 1 0 0;exp(lambda1*L) exp(lambda2*L) 0 0;0 0 1 1 ;0 0 exp(lambda1*L) exp(lambda2*L)];
bmat = [delta_T1_0;delta_T1_L;delta_T2_0;delta_T2_L];

G = Amat\bmat;
G1 = G(1);
G2 = G(2);
G3 = G(3);
G4 = G(4);

fL = (L*(G1*G4*lambda1 - G2*G3*lambda2) + (G2*G3 - G1*G4)*log((G1*exp(lambda1*L)+G2*exp(lambda2*L))/(G1+G2)))/(G1*G2*(lambda1-lambda2));
gL = (L*(G2*G3*lambda1 - G1*G4*lambda2) + (G1*G4 - G2*G3)*log((G3*exp(lambda1*L)+G4*exp(lambda2*L))/(G3+G4)))/(G3*G4*(lambda1-lambda2));

res(1) = log(delta_T1_L/delta_T1_0) - (U1*A_lm1*(1/C_h-1/C_c1) + U2*P2/C_h*fL);
res(2) = log(delta_T2_L/delta_T2_0) - (U2*A_lm2*(1/C_h-1/C_c2) + U1*P1/C_h*gL);

end

Apologies for posing a confusion.

Many thanks,

Chris

[Christiantyo Moeprodjo]

Thank you very much for your help Torsten! I really appreciate it.