Why doesn't Newton's method calculate the correct temperature

1 view (last 30 days)
Aryo Aryanapour
Aryo Aryanapour on 17 Jan 2022
Edited: Jon on 18 Jan 2022
function main
%UNTITLED Summary of this function goes here
% Detailed explanation goes here
x0 = 80;
eps = 0.4;
c_s = 5.67*10^(-8);
s1 = 0.250;
s2 = 0.015;
lamda1 = 0.35;
lamda2 = 22.7;
Tw_1 = 1200;
T_l = 10;
Tw_1 = 1200;
T_l = 10;
lambda_L = 25.12;
betha = 3.543*10^(-3);
v = 144.0*10^(-7);
L = 7;
B = 15;
A = B * L;
g = 9.81;
Pr = 0.7095;
Ra = @(x)Pr*(g*L^3*betha*(x-T_l))/v^2;
alpha_k = @(x)((0.825+(0.387*Ra(x)^1/6)/(1+(0.492/Pr)^9/16)^8/27)^2)*lambda_L/L;
func = @(x) (eps * c_s * x^4) + (alpha_k(x) + 1/(s1/lamda1+s2/lamda2)) * x - ((1/(s1/lamda1+s2/lamda2)) * Tw_1) + (alpha_k(x) * T_l);
xsol = newton_raphson(func,x0)
alpha_k = alpha_k(xsol)
alpha_ges = alpha_k + (eps *c_s*(xsol.^4-T_l.^4)/(xsol-T_l))
Q_ges = (1/((s1/lamda1+s2/lamda2)+1/alpha_ges))*A*(Tw_1-T_l)
(1/(s1/lamda1+s2/lamda2))*(Tw_1-xsol) % heat conduction
alpha_k * (xsol-T_l) + (eps*c_s*xsol^4) % convection + thermal radiation
end
function xsol = newton_raphson(func,x0)
%UNTITLED Summary of this function goes here
% Detailed explanation goes here
x(1) = x0;
maxiter = 200;
tol = 10^(-5);
for i = 1:maxiter
diff_xi = (func(x(i)+1e-6)-func(x(i)))*1e6;
if diff_xi < tol
fprintf('Pitfall hast occured a better initial guess\n');
return;
end
x(i+1) = x(i) - func(x(i))/diff_xi;
abs_error(i+1) = abs((x(i+1)-x(i))/x(i+1))*100;
if abs(x(i+1) - x(i)) < tol && (1/(s1/lamda1+s2/lamda2))*(Tw_1-x(i+1)) == alpha_k(x(i+1)) * (x(i+1)-T_l) + (eps*c_s*x(i+1)^4)
fprintf('The Root has converged at x = %.10f\n', x(i+1));
else
fprintf('Iteration no: %d,current guess x = %.10f, error = %.5f', i, x(i+1), abs_error(i+1));
end
end
xsol = x(end);
end
In the code above, the temperature, shown here as xsol, is to be calculated using the Newton method. Unfortunately, the heat transfer coefficient alpha_k also depends on the iterated xsol. I get a much too low temperature here. You can see that the code is not correct by the fact that heat conduction is not equal to heat convection+radiation.
(1/(s1/lamda1+s2/lamda2))*(Tw_1-xsol) % heat conduction
alpha_k * (xsol-T_l) + (eps*c_s*xsol^4) % convection + thermal radiation
Unfortunately, different values ​​come out here. I tried combining the if condition with an && to solve the problem. Unfortunately it does not work. I'm not very good at handling @functions at Matlab. I would be very grateful if someone could help me.
if abs(x(i+1) - x(i)) < tol && (1/(s1/lamda1+s2/lamda2))*(Tw_1-x(i+1)) == alpha_k(x(i+1)) * (x(i+1)-T_l) + (eps*c_s*x(i+1)^4)

Accepted Answer

Jon
Jon on 17 Jan 2022
Is there a reason you are using your own solver here rather than the MATLAB fzero (for function of one variable) or fsolve for systems of non-linear equations?
Unless there is a specific reason not to use those, I would try using the built in solver. If you still have convergence issues, then the problem is most likely with how you have formulated the equations, and not the solver. You can then focus on what the problem is there.
I just skimmed your code briefly, but one aspect that might be good to look at more closely, is whether in fact your heat transfer coefficient is getting updated using the current value of your solution. You have to look carefully at how the anonymous functions are defined and the variables are scoped to make sure it isn't actually remaining constant, just using the initial guess, x0
  9 Comments
Jon
Jon on 18 Jan 2022
The Newton-Raphson solves for a root (values that make function zero) of the function that you specify.
From your comments, it seems that you would like to find the temperature, x, that satisfies the energy balance:
conduction(x) = convection(x) + radiation(x)
or moving everything to the right hand side, you get the function you would like to find the root of
f(x) = -conduction(x) + convection(x) + radiation(x)
using the expressions for conduction, convection and radiation you have based upon the comments in your code you can define the function:
func = @(x) -(1/(s1/lamda1+s2/lamda2))*(Tw_1-x) + alpha_k(x) * (x-T_l) + (eps*c_s*x^4)
If I use this function, in place of yours, in your above code, I get:
Das Verfahren hat konvergiert bei x = 10.0000370032
Das Verfahren hat konvergiert bei x = 10.0000370032
xsol =
10.0000
alpha_k =
4.4982e+07
alpha_ges =
4.4982e+07
Q_ges =
1.7477e+05
ans =
1.6645e+03
ans =
1.6645e+03
So it seems like it converges ok and your checks are ok.
If you have reason to think the answer is not correct, you say you think it is too low, then you should probably recheck all of your expressions to make sure they are correct. Also, I would definitely be surprised if a celsius temperature should be used directly in an expression for the radiation. The radiation is typically in terms of Kelvin temperature to the 4th power, not Celsius as you have it.
Regarding your earlier termination criteria (original post)
if abs(x(i+1) - x(i)) < tol && (1/(s1/lamda1+s2/lamda2))*(Tw_1-x(i+1)) == alpha_k(x(i+1)) * (x(i+1)-T_l) + (eps*c_s*x(i+1)^4)
This condition will probably never be met as it requires that the energy balance be achieved exactly (you use ==) Due to round off and limited precision representation of floating point numbers you should never look for exact equality, instead you should test that the two values are within some tolerance of each other.
A few other points:
You use create a variable named eps in your program, this should be avoided as eps is already a built in function which gives the machine precision.
You don't need to create a function main that takes no arguments and does not return anything. You can just delete off the first line of code
function main
and the last end statement, and run it as a MATLAB script. Just type the name of the file on the command line to run the script.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!