34 views (last 30 days)

Show older comments

Basically, what the title asks. Code is here:

k = 180;

length = .05;

b = .01;

T0 = 180;

Tinf = 25;

h = 25;

Tsurr = 290;

epsilon = 0.9;

theta = atan((b/2)/length);

stef = 5.67*10^(-8);

xchange = b;

syms T1 T2 T3 T4 T5

eqn1 = (1-0.5*(xchange/length))*(T0 - T1) + (1-1.5*(xchange/length))*(T2-T1) + ((h*xchange^2)/(k*length*sin(theta)))*(Tinf - T1) + ((epsilon*stef*xchange^2)/(k*length*sin(theta)))*(Tsurr^4 - (T1+273)^4) == 0;

eqn2 = (1-1.5*(xchange/length))*(T1 - T2) + (1-2.5*(xchange/length))*(T3-T2) + ((h*xchange^2)/(k*length*sin(theta)))*(Tinf - T2) + ((epsilon*stef*xchange^2)/(k*length*sin(theta)))*(Tsurr^4 - (T2+273)^4) == 0;

eqn3 = (1-2.5*(xchange/length))*(T2 - T3) + (1-3.5*(xchange/length))*(T4-T3) + ((h*xchange^2)/(k*length*sin(theta)))*(Tinf - T3) + ((epsilon*stef*xchange^2)/(k*length*sin(theta)))*(Tsurr^4 - (T3+273)^4) == 0;

eqn4 = (1-3.5*(xchange/length))*(T3 - T4) + (1-4.5*(xchange/length))*(T5-T4) + ((h*xchange^2)/(k*length*sin(theta)))*(Tinf - T4) + ((epsilon*stef*xchange^2)/(k*length*sin(theta)))*(Tsurr^4 - (T4+273)^4) == 0;

eqn5 = 2*k*(xchange/2)*tan(theta)*((T4-T5)/xchange) + 2*h*((xchange/2)/cos(theta))*(Tinf - T5) + 2*epsilon*stef*((xchange/2)/cos(theta))*(Tsurr^4 - (T5 + 273)^4 ) == 0;

sol = solve([eqn1, eqn2, eqn3, eqn4, eqn5], [T1, T2, T3, T4, T5]);

T1sol = double(sol.T1)

T2sol = double(sol.T2)

T3sol = double(sol.T3)

T4sol = double(sol.T4)

T5sol = double(sol.T5)

Also, how long would this theoretically take to solve?

Thanks.

Matt J
on 26 Sep 2021 at 18:50

Well, one solution is clearly,

T1=T2=T3=T4=T5=(Tsurr-273)

as you can verify by direct subsitution.

Walter Roberson
on 27 Sep 2021 at 17:06

This is not one of the 1024 roots that my code comes up with.

Alan Weiss
on 27 Sep 2021 at 13:21

Edited: Alan Weiss
on 27 Sep 2021 at 13:23

Another point of view in case Matt's solution is not satisfactory: work completely in floating point computation rather than symbolic math computation. If you use the Problem-Based Optimization Workflow, you will basically have the same equations as before, simply using optimization variables rather than symbolic variables.

Good luck,

Alan Weiss

MATLAB mathematical toolbox documentation

Walter Roberson
on 27 Sep 2021 at 15:55

Edited: Walter Roberson
on 27 Sep 2021 at 17:04

k = 180;

length = .05;

b = .01;

T0 = 180;

Tinf = 25;

h = 25;

Tsurr = 290;

epsilon = 0.9;

theta = atan((b/2)/length);

stef = 5.67*10^(-8);

xchange = b;

syms T1 T2 T3 T4 T5

eqn1 = (1-0.5*(xchange/length))*(T0 - T1) + (1-1.5*(xchange/length))*(T2-T1) + ((h*xchange^2)/(k*length*sin(theta)))*(Tinf - T1) + ((epsilon*stef*xchange^2)/(k*length*sin(theta)))*(Tsurr^4 - (T1+273)^4) == 0;

eqn2 = (1-1.5*(xchange/length))*(T1 - T2) + (1-2.5*(xchange/length))*(T3-T2) + ((h*xchange^2)/(k*length*sin(theta)))*(Tinf - T2) + ((epsilon*stef*xchange^2)/(k*length*sin(theta)))*(Tsurr^4 - (T2+273)^4) == 0;

eqn3 = (1-2.5*(xchange/length))*(T2 - T3) + (1-3.5*(xchange/length))*(T4-T3) + ((h*xchange^2)/(k*length*sin(theta)))*(Tinf - T3) + ((epsilon*stef*xchange^2)/(k*length*sin(theta)))*(Tsurr^4 - (T3+273)^4) == 0;

eqn4 = (1-3.5*(xchange/length))*(T3 - T4) + (1-4.5*(xchange/length))*(T5-T4) + ((h*xchange^2)/(k*length*sin(theta)))*(Tinf - T4) + ((epsilon*stef*xchange^2)/(k*length*sin(theta)))*(Tsurr^4 - (T4+273)^4) == 0;

eqn5 = 2*k*(xchange/2)*tan(theta)*((T4-T5)/xchange) + 2*h*((xchange/2)/cos(theta))*(Tinf - T5) + 2*epsilon*stef*((xchange/2)/cos(theta))*(Tsurr^4 - (T5 + 273)^4 ) == 0;

E1 = [eqn1, eqn2, eqn3, eqn4, eqn5];

%{

partial2 = solve(E1(1), T2)

E2 = subs(E1(2:end), T2, partial2)

partial3 = solve(E2(1), T3)

E3 = subs(E2(2:end), T3, partial3)

partial4 = solve(E3(1), T4)

E4 = subs(E3(2:end), T4, partial4)

partial5 = solve(E4(1), T5)

E5 = subs(E4(2:end), T5, partial5)

%}

partials = solve(E1(1:4), [T2, T3, T4, T5]);

E5 = subs(lhs(E1(end))-rhs(E1(end)), partials);

At this point you can expand(E5) to get a plain polynomial in T1.

According to Maple, that plain polynomial is in degree 1024.

You should be able to take the following steps, but they take too much time for me to be able to run them and show the output here

tic

P = sym2poly(E5);

R = roots(P);

mask = imag(R) == 0;

T1_sol = R(mask);

T2_sol = double(subs(partials.T2, T1, T1_sol));

T3_sol = double(subs(partials.T3, T1, T1_sol));

T4_sol = double(subs(partials.T4, T1, T1_sol));

T5_sol = double(subs(partials.T5, T1, T1_sol));

toc

T1_sol

T2_sol

T3_sol

T4_sol

T5_sol

Walter Roberson
on 27 Sep 2021 at 17:05

T1_sol =

-9069.32334661296

-919.882898288149

177.024444908286

-24.4166503435627

T2_sol =

27738.4267958556

-2336.41842690512

174.081315968202

-287.461625512644

T3_sol =

7095829.16820471

-4126.23669567827

171.16801338191

-657.54975176523

T4_sol =

4.81617135399466e+16

-2960.81576251875

168.277608495836

-1280.43369676003

T5_sol =

3.06585398161644e+56

3425.69627315686

165.36383335676

-3127.23531356248

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

Start Hunting!