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

2 次查看（过去 30 天）
Ananthakrishnan
Ananthakrishnan 2024-5-5
评论： Ananthakrishnan 2024-5-6
Ran in:
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

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2024-5-5
Ran in:
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
  3 个评论
显示 1更早的评论
Torsten
Torsten 2024-5-5
编辑：Torsten 2024-5-5
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
Ananthakrishnan 2024-5-6
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

请先登录，再进行评论。

更多回答（1 个）

Venkat Siddarth Reddy
编辑：Venkat Siddarth Reddy 2024-5-5
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!

类别

在 Help Center 和 File Exchange 中查找有关 Electromagnetics 的更多信息

标签

产品


版本

R2024a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by