Code Running , Wrong answers . Information flow is wrong. Please help ! System of PDEs with coupled equations

Tried to solve a system of PDEs

Then Weq is depended on Td by equation :

where ϕ(i)=94.92211.*q(i)^4 - 76.5155.*q(i)^3 + 19.08577.*q(i)^2 + 0.03858.*q(1)+ 0.14278 ,ϕ(phi , relative humidity varying by input of q)

My code : its running but answers not coming as expected , probably cause information flow is not right. please show help me with it

clc
clear
close all

a = 0.05;%m
b = 0.025;%m
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 = 0.101; % kg/m
cd = 0.9; % 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].';
%4 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);

%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

plot(Td,Ta)


%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)=0.09;
 Ta(end)=297;
 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

Answer 1

[Torsten]

For Ta and w, you can only prescribe a value at z = 0, not at z = 1.

Further use upwind differencing for the first-order derivatives, not central differencing for Ta and w. Central differencing will introduce oscillations for the solutions.

The equations for Td and q are ordinary differential equations for which you are not allowed to prescribe any boundary values.