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

3 views (last 30 days)
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!
  2 Comments
Christiantyo Moeprodjo
Christiantyo Moeprodjo on 18 Nov 2022
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!

Sign in to comment.

Accepted Answer

Alan Stevens
Alan Stevens on 18 Nov 2022
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

More Answers (1)

Torsten
Torsten on 18 Nov 2022
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
  6 Comments
Christiantyo Moeprodjo
Christiantyo Moeprodjo on 29 Nov 2022
Edited: Torsten on 29 Nov 2022
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

Sign in to comment.

Categories

Find more on C Shared Library Integration in Help Center and File Exchange

Products


Release

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!