Hi everyone,
Im trying to solve a PDE using ode45. But the following error are ocurring: 'EDP must return a column vector'.
I put my code blow for you guys help me to find the error.
——————————————————————
global pipe Tin I Tg W T TW;T = [];Ts = [];TW = [];Tws = [];% condições iniciais
Iini = IC151(1,2);Tgini = TA075(1,2);Wini = FA061(1,2);Tini = TA072(1,2);T = Tini; pipe.Tref = 250; pipe.rho = 903 - 0.672 * pipe.Tref; pipe.Cp = 1820 + 3.76 * pipe.Tref; pipe.rhow = 7800; pipe.Cw = 450; pipe.ri = 0.009; pipe.ro = 0.013; pipe.di = pipe.ri * 2; pipe.do = pipe.ro * 2; pipe.Ai = pi * pipe.ri^2; pipe.Ao = pi * pipe.ro^2; pipe.ho = 11;pipe.hi=800; pipe.L=172; pipe.S=150; pipe.dX=pipe.L /pipe.S; pipe.n0=0.43; pipe.G=1.82; pipe.vs=Wini /pipe.Ai; pipe.gamma = (pipe.n0 * pipe.G) / (pipe.Ao * pipe.rhow * pipe.Cw);pipe.Tau1 = (pipe.Ai * pipe.rho * pipe.Cp) / (pi * pipe.di * pipe.hi);pipe.Tau12 = (pipe.Ao * pipe.rhow * pipe.Cw) / (pi * pipe.di * pipe.hi);pipe.Tau2 = (pipe.Ao * pipe.rhow * pipe.Cw) / (pi * pipe.do * pipe.ho);pipe.c = (pipe.vs * pipe.Tau1 * (1 + (pipe.Tau2 / pipe.Tau12)));TW = (Tini * pipe.Tau2) / (pipe.Tau12 + pipe.Tau2) + (Tgini * pipe.Tau12) / (pipe.Tau12 + pipe.Tau2) + (Iini * pipe.gamma * pipe.Tau12 * pipe.Tau2) / (pipe.Tau12 + pipe.Tau2);for(i=1:pipe.S) T(i+1) = pipe.gamma * pipe.Tau2 * Iini + Tgini + (T(1) - pipe.gamma * pipe.Tau2 * Iini - Tgini) * exp(-(i * pipe.dX) / pipe.c); TW(i+1) = (T(i+1) * pipe.Tau2) / (pipe.Tau12 + pipe.Tau2) + (Tgini * pipe.Tau12) / (pipe.Tau12 + pipe.Tau2) + (Iini * pipe.gamma * pipe.Tau12 * pipe.Tau2) / (pipe.Tau12 + pipe.Tau2);endTs = T(1,pipe.S+1);Tws = TW(1,pipe.S+1);n = 2;for k = 1:2481 I = IC151(n,2); % from workspace
Tin = TA072(n,2); % from workspace Tg = TA075(n,2); % from workspace W = FA061(n,2); % from workspace [dT dTW] = ode45(@EDP,[0 1],[T,TW],[]); n= n+1; end
——————————————————————
and the EDP file has:
——————————————————————-
function [dT, dTW] = EDP(t,x)global pipe W I Tg Tin T TW;T = x(1:pipe.S+1);T(1) = Tin;TW = x(pipe.S+1+1:2*(pipe.S+1));vs = W / pipe.Ai;dT(1) = 0;dTW(1) = 0;for(i=2:pipe.S+1) dT(i) = -vs * ((T(i) - T(i-1)) / pipe.dX) + ((TW(i) - T(i)) / pipe.Tau1); dTW(i) = pipe.gamma * I + ((Tg - TW(i)) / pipe.Tau2) - ((TW(i) - T(i)) / pipe.Tau12);end
——————————————————————
I would like to you help me.
best regards,
Gustavo
Best Answer