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

1 view (last 30 days)
Ananthakrishnan on 5 May 2024
Commented: Ananthakrishnan on 6 May 2024
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;
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
Array indices must be positive integers or logical values.

Error in solution>f (line 129)

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
dydt=zeros(length(y),1);
dqdt=zeros(Nz,1);
dwdt=zeros(Nz,1);
dTddt=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;
dTddz(i)=(Td(i+1)-Td(i-1))./2./dz;
dwdt(i)=-V.*dwdz(i)-(fd/2*A*rhoa).*dqdt(i);
end
end

Torsten on 5 May 2024
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;
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
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
dydt=zeros(length(y),1);
dqdt=zeros(Nz,1);
dwdt=zeros(Nz,1);
dTddt=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;
dTddz(i)=(Td(i+1)-Td(i-1))./2./dz;
dwdt(i)=-V.*dwdz(i)-(fd/2*A*rhoa).*dqdt(i);
end
end
3 CommentsShow 1 older commentHide 1 older comment
Torsten on 5 May 2024
Edited: Torsten on 5 May 2024
Fyi:
and
Ananthakrishnan on 6 May 2024
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

Venkat Siddarth Reddy on 5 May 2024
Edited: Venkat Siddarth Reddy on 5 May 2024
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!

Categories

Find more on Structural Mechanics in Help Center and File Exchange

R2024a

Community Treasure Hunt

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

Start Hunting!