I need HELP. Heat transfer ODE problem with bvp4c

I want to solve this equations using bcp4c but i can't find out what is wrong.

This is my code.

function dydx = odefun(x,y)

R = 2.5;ms =38.8 ;ma= 34.7;cpa= 1.007;cs=1.4;Uv=0.82;Ue=0.036291; Tamb=293;

dydx = zeros(2,1);

dydx(1) = (-1/(ms*cs))*((Uv*3.141569*(R^2)*(y(1)-y(2)))+(Ue*2*pi*R*(y(2)-Tamb)));

dydx(2) = (1/(ma*cpa))*(Uv*3.141569*(R^2)*(y(1)-y(2)));

end

bc = @(ya,yb) [ya(1) - 1073; yb(2) - 303]; %This are my boundary conditions ( Ts(x=0) = 1073 K ; Ta(x=40) = 303 K)

xmesh = linspace(0, 40, 50);

solinit = bvpinit(xmesh, [0 0]);

sol=bvp4c('DEdef2', bc, solinit);

figure;

plot(sol.x, sol.y);

[Jose Miguel Canales Huaman]

Please i really need help

Answer 1

[Alan Stevens]

You have an initial condition problem rather than a boundary condition problem (where you know values at the start and end x values), so bvp4c isn't appropriate. Use ode45 instead:

% Initial conditions
Ts0 = 1073;
Ta0 = 303;
y0 = [Ts0 Ta0];

xspan = [0 40];  % range of x values
% Call ode45
[x, y] = ode45(@odefun,xspan, y0);

% Plot results
plot(x,y(:,1),'r',x,y(:,2),'b'), grid
xlabel('x'), ylabel('y')

function dydx = odefun(~,y)
        R = 2.5;ms =38.8 ;ma= 34.7;cpa= 1.007;
        cs=1.4;Uv=0.82;Ue=0.036291; Tamb=293;
        dydx = [(-1/(ms*cs))*((Uv*3.141569*(R^2)*(y(1)-y(2)))+(Ue*2*pi*R*(y(2)-Tamb)));
               (1/(ma*cpa))*(Uv*3.141569*(R^2)*(y(1)-y(2)))];
end

[Jose Miguel Canales Huaman]

I thinks it is boundary value problem because the model describes a counter current heat exchanger where hot solid feed is one the left (X=0m) and cold air feed is on the right (x=40m).

I know the inlet conditions:

I want to know the profile tempetures of solid and air.

[Alan Stevens]

Apologies, I misread your conditions, thinking both were set at x = 0.

[Alan Stevens]

You seem to have the heat loss to ambient coming from the solid rather than the air (though this might not be the answer to your problem).

[Alan Stevens]

Thinking further, you need to solve for all of the temperatures simultaneously. The following does this (some explanation at the bottom):

Ts0 = 1073;
Ta0 = 303;
R = 2.5;ms =38.8 ;ma= 34.7;cpa= 1.007;
cs=1.4;Uv=0.82;Ue=0.036291; Tamb=293;

L = 40;
N = 40;
dx = L/N;
x = linspace(0, L, N);

k1 = Uv*pi*R^2*dx;
k2 = Ue*2*pi*R*dx;

M = zeros(2*N,2*N);
V = zeros(2*N,1);

M(1,1) = 1; M(2*N,2*N) = 1;
V(1) = Ts0; V(2*N) = Ta0;

for r = 2:N
    V(r) = 0;
    V(r+N-1) = k2/2*Tamb;
       c = r-1;
       M(r,c) = (k1/2 - ms*cs);
       M(r,c+1) = (k1/2 + ms*cs);
       M(r,c+N) = -k1/2;
       M(r,c+1+N) = -k1/2;
       M(r+N-1,c) = -k1/2;
       M(r+N-1,c+1) = -k1/2;
       M(r+N-1,c+N) = (k1/2 + k2/2 + ma*cpa);
       M(r+N-1,c+N+1) = (k1/2 +k2/2 - ma*cpa);
    
end

T = M\V;
Ts = T(1:N);
Ta = T(N+1:2*N);

plot(x,Ts,'r',x,Ta,'b'),grid
xlabel('x'), ylabel('Temperatures')
legend('Ts','Ta')

% Assume N temperature nodes with heat transfer in N-1 sections between.
%
% Solid:  ms*cs*(Ts(i+1) - Ts(i)) = -k1(Ts_av(i) - Ta_av(i))
% Air:    ma*cpa*(Ta(i)- Ta(i+1)) =  k1(Ts_av(i) - Ta_av(i)) - k2*(Ta_av(i) - Tamb)
%
% where Ts_av(i) = (Ts(i) + Ts(i+1))/2 and Ta_av(i) = (Ta(i) + Ta(i+1))/2
%
% The temperatures, Ts and Ta are represented by vector T of size 2Nx1
% where the first N are the values of Ts and the last N are the values of Ta.
%
% We assemble their coefficients in matrix M, and the other constants
% involving specified inlet temperatures and ambient temperature in 
% vector V.  Then we have  M*T = V, so we find T from T = M\V.

This results in:

NB I've switched the heat loss from solid_to_ambient to air_to_ambient as the former didn't seem to make much sense. However, you might have to switch that back if that's how it should really be.

[Jose Miguel Canales Huaman]

Its better than the results i got, but how can i switch that back?

Can you do it for me?

Here is the final graph that i'm looking for,:

Thanks for your help

[Alan Stevens]

I've tried that - it actually makes the comparison worse!

Have you checked the values of Uv and Ue? I don't have the information to do that.

[Jose Miguel Canales Huaman]

Can you give me your email, i can send you a pdf file. I'm working on a research and I want to take this model from a paper but i cant make it work.

Answer 2

[Alan Stevens]

Just upload the pdf as an attachment.

[Jose Miguel Canales Huaman]

Since page 95 to 119.

Thats the example i want to simulate,

[Alan Stevens]

I've now had a more detailed look at sections 6.2 and 6.3 of the dissertation. I've the following comments:

Page 100. Section 6.3.2

2nd paragraph says “Figure 6.2 illustrates the resulted temperature profiles of the solid and air along the cooler at two velocities of 0.2 and 0.5m/s.”

Figure 6.2 doesn’t illustrate temperature profiles, figure 6.3 does – a minor typographical error.

More seriously, using the data in table 6.1 the velocities (given by mdot/(rho*A) ) work out to be 1.37 m/s for air and 0.00089 m/s for the solid). Neither of these are in the range 0.2 to 0.5 m/s.

The paragraph goes on to say “Figure 6.2 (a) depicts the case of counter current flow while Fig. 6.2 (b) depicts the case of con-current flow.”

Allowing for the fact that 6.2 should be 6.3, there are no figures 6.3(a) and 6.3(b), just the one figure 6.3. There are tables 6.2(a) and 6.2(b), which are possibly what is meant, except that they don’t seem to refer to counter and con-current flows respectively, but to temperatures (a) and heat rates (b) for a counter-current case.

As figure 6.3 gives solid and air temperature traces over the length of the cooler, and table 6.3(a) gives the air temperature at x = 0, we can do an initial value solution going from x = 0 to x = L for both air and solid. The initial values being Ts(x=0) = 800°C, as before, with Ta(x=0) = 540°C. I’ve done this using equations 6.8 and 6.9 (with the sign of equation 6.9 changed to recognise that the air temperature will decrease going from x = 0 to x = L) and the data obtained from table 6.1, using MATLAB’s ode45 solver. The resulting temperature traces should match those of figure 6.3 if everything is consistent. They don’t! They are monstrously different – even going negative!

The above gives me no confidence that figure 6.3 represents the data provided in the tables. You are attempting the impossible by trying to reproduce it with that data!

[Jose Miguel Canales Huaman]

It seems to be a bad paper, with a lot of mistakes, but do you think the matematical model from equations 6.8 and 6.9 can be solve in matlab? , I also get negative values of tempeture, i was thinking that maybe I'm not writting a good code but now i have doubts about the mathematical model.

Another point is that both equations have to be solved simultaneously, because variable Ta and Ts change as x change, i dont think ode45 can solve it. Thanks for your help. Do you know a blog or website where i can find more help?

O maybe do you know a person that can solve equiations 6.8 and 6.9 in mathlab without negative results?

Thanks a lot :D

[Alan Stevens]

Here's a simple Euler code that solves the equations 6.8 and 6.9, using your original data values.

Ts0 = 800;
Ta0 = 30;
R = 2.5;ms =38.8889 ;ma= 34.7;cpa= 1.007;
cs=1.4;Uv=0.82;Ue=0.036291; Tamb=20;

L = 40;
N = 18;
dx = L/N;
x = linspace(0, L, N)/L;

k1 = Uv*pi*R^2*dx;
k2 = Ue*2*pi*R*dx;

M = zeros(2*N,2*N);
V = zeros(2*N,1);

M(1,1) = 1; M(2*N,2*N) = 1;
V(1) = Ts0; V(2*N) = Ta0;

for r = 2:N
    V(r) = k2*Tamb;
    V(r+N-1) = 0;
       c = r-1;
       M(r,c) = (k1 - ms*cs);
       M(r,c+1) = ms*cs;
       M(r,c+N) = -k1+k2;
       M(r+N-1,c) = -k1;
       M(r+N-1,c+N) = (k1 + ma*cpa);
       M(r+N-1,c+N+1) = - ma*cpa;
    
end

T = M\V;
Ts = T(1:N);
Ta = T(N+1:2*N);

plot(x,Ts,'rs-',x,Ta,'bo-'),grid
axis([0 1 0 800])
xlabel('x'), ylabel('Temperatures')
legend('Ts','Ta')

% Assume N temperature nodes with heat transfer in N-1 sections between.
%
% Solid:  ms*cs*(Ts(i+1) - Ts(i)) = -k1(Ts(i) - Ta(i)) - k2*(Ta(i) - Tamb)
% Air:    ma*cpa*(Ta(i)- Ta(i+1)) =  k1(Ts(i) - Ta(i)) 
%
%
% The temperatures, Ts and Ta are represented by vector T of size 2Nx1
% where the first N are the values of Ts and the last N are the values of Ta.
%
% We assemble their coefficients in matrix M, and the other constants
% involving specified inlet temperatures and ambient temperature in 
% vector V.  Then we have  M*T = V, so we find T from T = M\V.

% Ts(1) = Ts0; Ta(2*N) = Ta0;

[Jose Miguel Canales Huaman]

Why it is not possible to solve this equations as a bvp?

and another question, How can i use your code but without the change of minus that you mention before.

[Alan Stevens]

You could try to solve it using bvp4c. I never use it, finding its syntax cumbersome.

The latest code I uploaded doesn't have any sign change - that was in the code (not uploaded) where I integrated from x=0 to x=L for both equations. The last code uploaded basically explicitly discretises the length into sections, uses a simple Euler type solution for each section, assembles them in a matrix and solves all sections simultaneously using the matrix division.

[Jose Miguel Canales Huaman]

Thanks a lot, I`m going to test your results with experimental data, thanks :D

[Jose Miguel Canales Huaman]

% Air: ma*cpa*(Ta(i)- Ta(i+1)) = k1(Ts(i) - Ta(i))

In this part, why did you type Ta(1) - Ta(i+1), because diference of tempetures is T final - T inicial, so it should be Ta(i+1) - Ta(1). Or am I misunderstanding?

[Alan Stevens]

Think of the physics of the situation. Look at the right-hand side of the air energy equation. It is positive, because k1 is positive and Ts(i) is greater than Ta(i). This means the left hand side will be positive. Now we know that the air temperature will be higher towards the left because the air picks up heat as it flows from right to left. Ta(i) is further to the left than Ta(i+1), so the left hand side must have Ta(i) - Ta(i+1), not Ta(i+1) - Ta(i). You can test this by changing the sign of the two values of ma*cpa in my last program. You will get nonsensical results.

That said, I still think the equations are not strictly correct because the heat loss to atmosphere is taken off the solid energy equation. I think it should come off the air energy equation instead. Numerically, this doesn't make a great difference with your data though, because the effective heat transfer coefficient to atmospere is relatively very small.