1 view (last 30 days)

Show older comments

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)) ;

end

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]

end

double(T_p_1)

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

syms(symvar(eqn));

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
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

Konrad
on 12 Aug 2021

Hi,

One quick guess:

in your equation (eqn(number_of_nodes) = ...) you are raising to the power of .24, which means you are taking the 4.166..th 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

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

Start Hunting!