Azzera filtri
Azzera filtri

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

3 visualizzazioni (ultimi 30 giorni)
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

Risposte (1)

Torsten
Torsten il 6 Mag 2024
Spostato: Torsten il 6 Mag 2024
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.

Prodotti


Release

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by