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

3 ビュー (過去 30 日間)
Aryo Aryanapour
Aryo Aryanapour 2022 年 1 月 17 日
編集済み: Jon 2022 年 1 月 18 日
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)

採用された回答

Jon
Jon 2022 年 1 月 17 日
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 件のコメント
Jon
Jon 2022 年 1 月 17 日
It is hard to follow what your code is trying to do. I'm sorry I don't have time now to go through it and try to understand each step. Maybe someone else can assist you further.
Jon
Jon 2022 年 1 月 18 日
編集済み: Jon 2022 年 1 月 18 日
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.

サインインしてコメントする。

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeProgramming についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by