How do you solve a System of 11 Non-Linear Equations Inter-linked to one another?
3 views (last 30 days)
Show older comments
Christiantyo Moeprodjo
on 17 Nov 2022
Commented: Christiantyo Moeprodjo
on 4 Dec 2022
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
Accepted Answer
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])
disp(fval)
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
0 Comments
More Answers (1)
Torsten
on 18 Nov 2022
U0 = [663 255];
options = optimset('TolX',1e-12,'TolFun',1e-12);
format long
U = fsolve(@fun,U0,options)
fun(U)
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
See Also
Categories
Find more on C Shared Library Integration in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!