could someone please tell me why my code returns complex answers while it is about room temperature ?

1 view (last 30 days)
clc , clear
% inputs :
ho = 34 ;
k1 = 25;
cp1 = 10 ;
c1 = 12682.7 ;
ro = 3410.9 ;
q = 0 ;
desired_time=25; %sec
dt = 0.5 ;
n = desired_time/dt ;
wall_thickness = 219; %mm
number_of_nodes = 5 ;
dx = ( wall_thickness / (number_of_nodes - 1) )*10^-3 ;
outside_temperature = 10 ;
To = sym(outside_temperature) ;
initial_room_temperature = 25 ;
initial_wall_temperature = 20 ;
% equations :
T_p_1_n = zeros(1,number_of_nodes+1)+ initial_wall_temperature ;
T_p_1_n(1,number_of_nodes+1) = initial_room_temperature ;
T_p_1 = sym(T_p_1_n);
T_p_a = sym(zeros(n,number_of_nodes+1)) ;
T_p = sym('T_p%d' ,[1,number_of_nodes+1],'real');
eqn = zeros(1,number_of_nodes+1) ;
eqn = sym(eqn) ;
for o=1:n
eqn(1) = ho * (To - T_p(1)) - (k1 / dx)*(T_p(1) - T_p(2)) - ((ro * cp1 * (dx/2))/dt) * (T_p(1) -T_p_1(1)) ; % outside of the wall node
for i = 2:number_of_nodes-1 %inside wall nodes
eqn(i) = (k1/dx) * (T_p(i-1) - T_p(i)) - (k1/dx)*(T_p(i) - T_p(i+1)) - ((ro * cp1 * dx )/dt) * (T_p(i)-T_p_1(i)) ;
eqn(number_of_nodes) = (k1/dx) / ( T_p(number_of_nodes-1) - T_p(number_of_nodes)) - (2.3* ((T_p(number_of_nodes) - T_p(number_of_nodes+1))^0.24)) * (T_p(number_of_nodes) - T_p(number_of_nodes+1)) - ((ro * cp1 * (dx/2))/dt)*(T_p(number_of_nodes)-T_p_1(number_of_nodes)) ; % inside of the wall node
eqn(number_of_nodes+1) = q - (28.9805) * (2.3* ((T_p(number_of_nodes) - T_p(number_of_nodes+1))^0.24)) * (T_p(number_of_nodes+1)-T_p(number_of_nodes)) - (c1/dt) * (T_p(number_of_nodes+1)-T_p_1(number_of_nodes+1)) ; % room air equation
[x1 , x2 , x3 , x4 , x5 , x6] = vpasolve( eqn,symvar(eqn) )
T_p_1 = [x1 , x2 , x3 , x4 , x5 , x6]

Accepted Answer

Walter Roberson
Walter Roberson on 12 Aug 2021
[x1 , x2 , x3 , x4 , x5 , x6] = vpasolve( eqn,symvar(eqn) )
if you breakpoint at that line, and run to there, then you can do
sol1_4 = solve(eqn(1:4), [T_p1, T_p2, T_p3, T_p4]) %easy
neqn = simplify(subs(eqn(5:end), sol1_4))
It is then possible to solve neqn(1) for T_p5 -- or at least Maple can do it. The result is T_p5 expressed in terms of the roots of a polynomial of degree 56 with non-zero coefficients at degree 56, 50, 31, 25, and 0.
Now you could hope for a moment to factor into a polynomial of degree 31 times a polynomial of degree 25, but that would get you entries with degree 56, 31, 25, and 0... but no entries of degree 50. And if you try to include (polynomial of degree 25)^2 then you will start to get additional powers; conjugate roots can get you degree 50 but not the 25 needed to add to 31 to get 56. Anyhow, you probably cannot factor.
Once you have the root() form of the polynomial of degree 56, you can substitute it into neqn(2) to get the equation that needs to hold for T_p6 . You are not going to be able to solve that symbolically. If you try to solve it numerically... you have problems. If you plot the imaginary part, it looks like there might be only one place the imaginary part is 0, near 20.749... but at that location, the real part is over 100000.
It appears that all of the roots for T_p6 are imaginary.
Walter Roberson
Walter Roberson on 12 Aug 2021
Edited: Walter Roberson on 12 Aug 2021
Now, these plots do depend upon MATLAB having solved equations correctly... namely the linear equations in the first 4 parts. But even if it got the exact values of the linear solutions wrong, it would have to have gotten the form of the solutions wrong for it to have ruined the possibility of a solution for the remaining two equations. And my work with the equations in Maple shows that MATLAB did not get the form of the solutions wrong for the first four equations.
The surface plots do not depend upon MATLAB being able to solve the remaining two equations: only upon it being able to evaluate the equations at points and graph the results correctly. MATLAB is not without its flaws in graphing equations that require high precision because they have terms that nearly balance, but this is not one of those situations.
At this point, we are faced with one of the following possibilities:
  • the laws of thermodynamics might be wrong
  • fundamental Mathematics might be wrong
  • you might have made a mistake in deriving the equations
  • you might have made a mistake in transcribing the equations into MATLAB
  • MATLAB and Maple might both have major major mistakes in calculating the form of solutions to simple linear equations
  • MATLAB and Maple might both have significant problems in evaluating the objective function obj() that is calculated here

Sign in to comment.

More Answers (1)

Konrad on 12 Aug 2021
One quick guess:
in your equation (eqn(number_of_nodes) = ...) you are raising to the power of .24, which means you are taking the root. If the term inside the parentheses
(T_p(number_of_nodes) - T_p(number_of_nodes+1))^.24
is negative, you get a complex result.
Best, Konrad
roham farhadi
roham farhadi on 12 Aug 2021
sorry for bad explanations ,
this code use explicit finite-difference method to find transient room air and walls temperature so the basic is in each time step a system of equations must be solved to find T_ P , and then T_P_1 (which is initial conditions ) must be replaced with T_p for the next time step , and this continued to reach the desired time (like 25 sec ) , this is what the main loop do.

Sign in to comment.


Find more on Thermal Analysis 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!

Translated by