how can i code energy balance eq in matlab

i need to model these two eqns of heat exchanger

dT_Co/dt = ((M_C*(T_Ci-T_Co))/(rho_C*V_C))+(UA_i*(T_Ho-T_Co)/(rho_C*V_C*C_Pmc))

dT_Ho/dt = ((M_H*(T_Hi-T_Ho))/(rho_H*V_H))+(UA_i*(T_Co-T_Ho)/(rho_H*V_H*C_Pmh))

where all the inputs are known to me except T_Co and T_Ho in both the equations. the resultant T_Co and T_Ho will be the inputs for finding the differential eqns. can some one please help me with this. Thanks in advance.

[John D'Errico]

help ode45

help dsolve

Answer 1

[Stephan]

There is an analytical solution found by dsolve:

syms T_Co(t) T_Ho(t) M_C M_H T_Ci T_Hi rho_C rho_H V_C V_H UA_i C_Pmc C_Pmh
eq(1) = diff(T_Co,t) == ((M_C*(T_Ci-T_Co))/(rho_C*V_C))+(UA_i*(T_Ho-T_Co)/(rho_C*V_C*C_Pmc));
eq(2) = diff(T_Ho,t) == ((M_H*(T_Hi-T_Ho))/(rho_H*V_H))+(UA_i*(T_Co-T_Ho)/(rho_H*V_H*C_Pmh));
sol = dsolve(eq)
sol_T_Ho = sol.T_Ho
sol_T_Co = sol.T_Co

Read the link to see how you can use initial conditions with this.

If you want to explore this system numerical use:

syms T_Co(t) T_Ho(t) M_C M_H T_Ci T_Hi rho_C rho_H V_C V_H UA_i C_Pmc C_Pmh
eq(1) = diff(T_Co,t) == ((M_C*(T_Ci-T_Co))/(rho_C*V_C))+(UA_i*(T_Ho-T_Co)/(rho_C*V_C*C_Pmc));
eq(2) = diff(T_Ho,t) == ((M_H*(T_Hi-T_Ho))/(rho_H*V_H))+(UA_i*(T_Co-T_Ho)/(rho_H*V_H*C_Pmh));

% create a function handle from the symbolic system
[eqs,vars] = odeToVectorField(eq);
fun = matlabFunction(eqs,'vars',{'t','Y','M_C', 'M_H', 'T_Ci', 'T_Hi', 'rho_C', 'rho_H',...
    'V_C', 'V_H', 'UA_i', 'C_Pmc', 'C_Pmh'})

% Making some fantasy values for numeric calculation
M_C_num = 3;
M_H_num = 2;
T_Ci_num = 1;
T_Hi_num = 7;
rho_C_num = 2;
rho_H_num = 9;
V_C_num = 2;
V_H_num = 5;
UA_i_num = 4;
C_Pmc_num = 3;
C_Pmh_num = 2;
y0 = [100 0];
tspan = [0 200];

% solve the system
[t,res] = ode45(@(t,Y)fun(t,Y,M_C_num,M_H_num,T_Ci_num,T_Hi_num,rho_C_num,rho_H_num,V_C_num,...
    V_H_num,UA_i_num,C_Pmc_num,C_Pmh_num),tspan,y0);

% plot results
figure(1)
plot(t,res)
figure(2)
plot3(res(:,1),res(:,2),t,'-r','LineWidth',2)
grid on
xlabel('THo')
ylabel('TCo')
zlabel('Time')

which results in:

and:

[vishnuvardhan naidu tanga]

hello Stephan, thank you for the solution. i have already tried this but there is a license issue with my software. I will look into it again. thank you.

with regards,vishnuvardhan

[Stephan]

I suspect that you do not have symbolic toolbox. In this case:

fun = @(t,Y,M_C,M_H,T_Ci,T_Hi,rho_C,rho_H,V_C,V_H,UA_i,C_Pmc,C_Pmh)...
    [(M_H.*(T_Hi-Y(1)))./(V_H.*rho_H)-(UA_i.*(Y(1)-Y(2)))./(C_Pmh.*V_H.*rho_H);...
    (M_C.*(T_Ci-Y(2)))./(V_C.*rho_C)+(UA_i.*(Y(1)-Y(2)))./(C_Pmc.*V_C.*rho_C)];

M_C_num = 3;
M_H_num = 2;
T_Ci_num = 1;
T_Hi_num = 7;
rho_C_num = 2;
rho_H_num = 9;
V_C_num = 2;
V_H_num = 5;
UA_i_num = 4;
C_Pmc_num = 3;
C_Pmh_num = 2;
y0 = [100 0];
tspan = [0; 200];

[t,res] = ode45(@(t,Y)fun(t,Y,M_C_num,M_H_num,T_Ci_num,T_Hi_num,rho_C_num,rho_H_num,V_C_num,...
    V_H_num,UA_i_num,C_Pmc_num,C_Pmh_num),tspan,y0);

figure(1)
plot(t,res)

figure(2)
plot3(res(:,1),res(:,2),t,'-r','LineWidth',2)
grid on
xlabel('THo')
ylabel('TCo')
zlabel('Time')

[vishnuvardhan naidu tanga]

hello Stephan, thank you for the solution. the code now runs in my system.

[vishnuvardhan naidu tanga]

hello Stephan, if i want to use this code as a function how can i do it? i am trying to implement this as a function in simulink, but it keeps on getting an error message as the deriavative t is not defined.