Trying to solve system of PDE by spacial discretization. Errors are coming frequent . what am i doing wrong ?

clc
clear
close all

a = 10;
b = 10;
Lp = 2*b + 2*(sqrt(b^2 + (a*pi)^2))*((3 + (2*b/a*pi)^2)/(4 + (2*b/a*pi)^2));
A = 2*a*b;
Pa = 101325; % Pa
Ta = 30; % Celsius
D = (2.256/Pa)*((Ta+273.13)/256)^1.81;
V = 2;
rhoa = 1.15; % kg/m^3
Sha = 2.45;
Ky = rhoa*Sha*D*Lp/4*A;
Nua = 2.45;
ka = 0.0321; % W/(m K)
h = Nua*ka*Lp/4*A;
had = 0.9; % J/kg K (specific heat of adsorption)
fd = 0.914; % kg/m
fm = 12; % kg/m
cd = 0.85; % J/kg K
cw = 4.19; % J/kg K
cm = 0.9; % J/kg K
ca = 1.005; % J/kg K
cv = 2.028; % J/kg K


w10=0.0133 % intital values 

w10 = 0.0133

Ta10=303

Ta10 = 303

Td10=303

Td10 = 303

RH10=0.5 

RH10 = 0.5000

q10=0.01

q10 = 0.0100

%increment for lenght
Nz=100;
z=linspace(0,1,Nz);
dz=z(2)-z(1);


%time
t=0:3:299;%secs

%inital conditions all zeros
ICi=zeros(1,Nz);
IC=[ICi,ICi,ICi,ICi];
%6 variables os IC must be 4 c1-c6

%odesolver
[t,y]=ode15s(@f,t,IC,[],Lp,A,Pa,Ta,V,rhoa,Ky,ka,h,had,fd,fm,cd,cw,cm,ca,cv,w10,Ta10,Td10,q10,Nz,dz);

Array indices must be positive integers or logical values.

Error in solution>f (line 129)
dTadt(i)=-V.*dTadz(i)-((fd(cd+q(i).*cw)+fm*cm)/(A*rhoa*(ca+w(i)*cv))).*dTddt(i)+(ky*P*had)/(A*rhoa*(ca+w*cv))*(w(i)-weq(i));

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode15s (line 148)
odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

%re calculation
%define value

q=y(:,1:Nz);
w=y(:,Nz+1:2*Nz);
Td=y(:,2*Nz+1:3*Nz);
Ta=y(:,3*Nz+1:4*Nz);


%bounday cond
q(:,1)=q10 %kg/kgda
w(:,1)=w10; %celcius
Td(:,1)=Td10; %celcius
Ta(:,1)=Ta10

q(:,end)=w(:,end-1);
w(:,end)=w(:,end-1);
Ta(:,end)=Ta(:,end-1);
Td(:,end)=Td(:,end-1);

%plot


%PDE-ODe 
function dydt=f(~,y,Lp,A,Pa,Ta,V,rhoa,Ky,ka,h,had,fd,fm,cd,cw,cm,ca,cv,w10,Ta10,Td10,q10,Nz,dz)
 
dydt=zeros(length(y),1);

 dqdt=zeros(Nz,1);
 dwdt=zeros(Nz,1);
 dTddt=zeros(Nz,1);
 dTadt=zeros(Nz,1);


 %define value
 q=y(1:Nz);
 w=y(Nz+1:2*Nz);
 Td=y(2*Nz+1:3*Nz);
 Ta=y(3*Nz+1:4*Nz);

 %bounday condition

 q(1)=q10;%kg/kg
 w(1)=w10;%30 celcius
 Td(1)=Td10;%30celcius
 Ta(1)=Ta10; 

 q(end)=q(end-1);
 w(end)=w(end-1);
 Ta(end)=Ta(end-1);
 Td(end)=Td(end-1); 
 

 %interior
 for i=2:Nz-1
    pws(i)=exp(23.196-(3816.44/(Td(i)-46.13)));

    
    RH(i)=94.92211.*q(i)^4 - 76.5155.*q(i)^3 + 19.08577.*q(i)^2 + 0.03858.*q(1)+ 0.14278;
    weq(i)=0.62188*RH(i)./((Pa./pws(i))-RH(i));
    dqdt(i)=-2*Ky*Pa/fd*(weq(i)-w(i)); % one of the  sysetm equations but its an ODE 



    dqdz(i)=(q(i+1)-q(i-1))./2./dz;

    dwdz(i)=(w(i+1)-w(i-1))./2./dz;
    dTadz(i)=(Ta(i+1)-Ta(i-1))./2./dz;
    dTddz(i)=(Td(i+1)-Td(i-1))./2./dz;

    dwdt(i)=-V.*dwdz(i)-(fd/2*A*rhoa).*dqdt(i);
    dTddt(i)=-((2*h*Lp)/(fd*(cd+q(i)*cw)+fm*cm)).*(Td(i)-Ta(i))-((2*Ky*Lp*had)/(fd*(cd+q(i)*cw)+fm*cm)).*(weq(i)-w(i))-((2*Ky*Lp*cv)/(fd*(cd+q(i)*cw)+fm*cm)).*(weq(i)-w(i));
    dTadt(i)=-V.*dTadz(i)-((fd(cd+q(i).*cw)+fm*cm)/(A*rhoa*(ca+w(i)*cv))).*dTddt(i)+(ky*P*had)/(A*rhoa*(ca+w*cv))*(w(i)-weq(i));
    
 end
 dydt=[dqdt,dwdt,dTddt,dTadt]
 end

Answer 1

[Torsten]

Technically, your code works now. But you get NaN values for some dydt's. The reason might be that you start with a complete 0 vector of initial conditions and thus divide by 0 somwehere. But it's time you start debugging ...

clc
clear
close all

a = 10;
b = 10;
Lp = 2*b + 2*(sqrt(b^2 + (a*pi)^2))*((3 + (2*b/a*pi)^2)/(4 + (2*b/a*pi)^2));
A = 2*a*b;
Pa = 101325; % Pa
Ta = 30; % Celsius
D = (2.256/Pa)*((Ta+273.13)/256)^1.81;
V = 2;
rhoa = 1.15; % kg/m^3
Sha = 2.45;
Ky = rhoa*Sha*D*Lp/4*A;
Nua = 2.45;
ka = 0.0321; % W/(m K)
h = Nua*ka*Lp/4*A;
had = 0.9; % J/kg K (specific heat of adsorption)
fd = 0.914; % kg/m
fm = 12; % kg/m
cd = 0.85; % J/kg K
cw = 4.19; % J/kg K
cm = 0.9; % J/kg K
ca = 1.005; % J/kg K
cv = 2.028; % J/kg K


w10=0.0133; % intital values 
Ta10=303;
Td10=303;
RH10=0.5 ;
q10=0.01;




%increment for lenght
Nz=100;
z=linspace(0,1,Nz);
dz=z(2)-z(1);


%time
t=0:3:299;%secs

%inital conditions all zeros
ICi=zeros(1,Nz);
IC=[ICi,ICi,ICi,ICi].';
%6 variables os IC must be 4 c1-c6

%odesolver
[t,y]=ode15s(@(t,y)f(t,y,Lp,A,Pa,Ta,V,rhoa,Ky,ka,h,had,fd,fm,cd,cw,cm,ca,cv,w10,Ta10,Td10,q10,Nz,dz),t,IC);

Warning: Failure at t=1.904303e-07. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (4.235165e-22) at time t.

%re calculation
%define value

q=y(:,1:Nz);
w=y(:,Nz+1:2*Nz);
Td=y(:,2*Nz+1:3*Nz);
Ta=y(:,3*Nz+1:4*Nz);


%bounday cond
q(:,1)=q10; %kg/kgda
w(:,1)=w10; %celcius
Td(:,1)=Td10; %celcius
Ta(:,1)=Ta10;

q(:,end)=w(:,end-1);
w(:,end)=w(:,end-1);
Ta(:,end)=Ta(:,end-1);
Td(:,end)=Td(:,end-1);

%plot


%PDE-ODe 
function dydt=f(~,y,Lp,A,Pa,Ta,V,rhoa,Ky,ka,h,had,fd,fm,cd,cw,cm,ca,cv,w10,Ta10,Td10,q10,Nz,dz)
 
dydt=zeros(length(y),1);

 dqdt=zeros(Nz,1);
 dwdt=zeros(Nz,1);
 dTddt=zeros(Nz,1);
 dTadt=zeros(Nz,1);


 %define value
 q=y(1:Nz);
 w=y(Nz+1:2*Nz);
 Td=y(2*Nz+1:3*Nz);
 Ta=y(3*Nz+1:4*Nz);

 %bounday condition

 q(1)=q10;%kg/kg
 w(1)=w10;%30 celcius
 Td(1)=Td10;%30celcius
 Ta(1)=Ta10; 

 q(end)=q(end-1);
 w(end)=w(end-1);
 Ta(end)=Ta(end-1);
 Td(end)=Td(end-1); 
 

 %interior
 for i=2:Nz-1
    pws(i)=exp(23.196-(3816.44/(Td(i)-46.13)));

    
    RH(i)=94.92211.*q(i)^4 - 76.5155.*q(i)^3 + 19.08577.*q(i)^2 + 0.03858.*q(1)+ 0.14278;
    weq(i)=0.62188*RH(i)./((Pa./pws(i))-RH(i));
    dqdt(i)=-2*Ky*Pa/fd*(weq(i)-w(i)); % one of the  sysetm equations but its an ODE 



    dqdz(i)=(q(i+1)-q(i-1))./2./dz;

    dwdz(i)=(w(i+1)-w(i-1))./2./dz;
    dTadz(i)=(Ta(i+1)-Ta(i-1))./2./dz;
    dTddz(i)=(Td(i+1)-Td(i-1))./2./dz;

    dwdt(i)=-V.*dwdz(i)-(fd/2*A*rhoa).*dqdt(i);
    dTddt(i)=-((2*h*Lp)/(fd*(cd+q(i)*cw)+fm*cm)).*(Td(i)-Ta(i))-((2*Ky*Lp*had)/(fd*(cd+q(i)*cw)+fm*cm)).*(weq(i)-w(i))-((2*Ky*Lp*cv)/(fd*(cd+q(i)*cw)+fm*cm)).*(weq(i)-w(i));
    dTadt(i)=-V.*dTadz(i)-((fd*(cd+q(i).*cw)+fm*cm)/(A*rhoa*(ca+w(i)*cv))).*dTddt(i)+(Ky*Pa*had)/(A*rhoa*(ca+w(i)*cv))*(w(i)-weq(i));
    
 end
 dydt=[dqdt;dwdt;dTddt;dTadt];
 end

[Ananthakrishnan]

Thanks a lot, i need to figure out dimensions in my constants rn.

[Torsten]

Fyi:

The changes I made were

dTadt(i)=-V.*dTadz(i)-((fd*(cd+q(i).*cw)+fm*cm)/(A*rhoa*(ca+w(i)*cv))).*dTddt(i)+(Ky*Pa*had)/(A*rhoa*(ca+w(i)*cv))*(w(i)-weq(i));

instead of

dTadt(i)=-V.*dTadz(i)-((fd(cd+q(i).*cw)+fm*cm)/(A*rhoa*(ca+w(i)*cv))).*dTddt(i)+(ky*P*had)/(A*rhoa*(ca+w*cv))*(w(i)-weq(i));

and

dydt=[dqdt;dwdt;dTddt;dTadt];

instead of

dydt=[dqdt,dwdt,dTddt,dTadt] 

[Ananthakrishnan]

Hello , i have done the corrections and updated the code and it got working, but the information flow my eqquations seemed to correct. i have posted in the forum regarding this . if possible can you help me with that

Answer 2

[Venkat Siddarth Reddy]

Hi Ananthakrishnan,

I think you are missing an arithmetic operation in the last line of for loop in function f next to the variable fd i.e.

Since fd is a scalar value and the value of equation in the parenthesis next to it is a fractional. It throws an error.

I hope it helps!