finite difference implicit method
显示 更早的评论
Hi everyone, I have written this code but I do not know why Matlab does not read the if condition. It suppose to use different variable for (alfa) when it is reach N= 33, 66.
Could you please help me with that.
clear all
clc
ho=50;
hi=10;
N=99; %Number of space intervals
M=2000; %number of time steps
Th=196; %thickness of wall in meter;;
t_0=0; %initial time
t_f=3600; %final time (sec)
To=34; %exterior wall temperature
k_b=0.7; %thermal conductivity of bricks
roh_b=1600; %density of bricks (kg/m^3)
Cp_b=840; %Cp of bricks (J/kg-K)
T_in=22; %Inner wall temperature
abs=0.7;
dx = (Th/(N-1))*0.001; %size of the mesh
% dx=0.002;
x = 0:dx:N-1;
x = x'; %inverse of mesh
bi=hi*dx/k_b; %coefficient of terms
bo=ho*dx/k_b;
%n-Octadecane PCM
k_ps=0.358;
k_pl=0.148;
Cp_ps=1.934;
Cp_pl=2.196;
roh_ps=865;
roh_pl=780;
lamda=243.5;
dT=1;
Tm=25;
Gt=solarradiationcode;
%initial condition
u(:,1)=ones(N,1)*15;
for index=1:length(Gt)
% index=1;
for j=1:M %time step
V_Te=To*(abs*Gt/ho);
Te=V_Te(index);
dt =3;
t = t_0:dt:t_f;
alfa= k_b/(Cp_b*roh_b);
% define the ratio r
r=alfa*dt/dx^2; %Fourier number
Beta=[2*r*bo*Te; zeros(N-2,1); 2*r*bi*T_in];
A = zeros(N,N);
A(1,1) = 1+2*r+2*r*bo;
A(1,2) = -2*r;
for i=2:N-1
if N/3> i || i>2*N/3
Cp=Cp_b;
roh=roh_b;
k=k_b;
alfa= k/(Cp*roh);
r=alfa*dt/dx^2;
end
if N/3==i && i==2*N/3
Cp=2*Cp_b*Cp_pl/(Cp_pl+Cp_b);
roh=2*roh_pl*roh_b/(roh_pl+roh_b);
k=2*k_b*k_pl/(k_pl+k_b);
alfa= k/(Cp*roh);
r=alfa*dt/dx^2;
end
if N/3<i && i<2*N/3
T=u(i,j);
if T<Tm
Cp=Cp_ps;
roh=roh_ps;
k=k_ps;
elseif Tm<T && T<Tm+dT
Cp= ((Cp_ps+Cp_pl)/2)+((roh_ps+roh_pl)/2)*(lamda/dT);
elseif T>Tm+dT
Cp=Cp_pl;
roh=roh_pl;
k=k_pl;
end
alfa=k/(Cp*roh);
r=alfa*dt/dx^2;
end
A(i,i-1) = -r;
A(i,i) = 1+2*r;
A(i,i+1) = -r;
end
alfa= k_b/(Cp_b*roh_b);
% define the ratio r
r=alfa*dt/dx^2; %Fourier number
%A(N,N+1) = 0;
A(N,N-1) = -2*r;
A(N,N) = 1+2*r+2*r*bi;
[L,U] = tridia(A);
[v] = lower_solver(L,u(:,j)+Beta);
u(:,j+1)= uper_solve(U,v);
end
end
plot(u(:,j+1))
回答（1 个）
the cyclist
2014-5-1
编辑：the cyclist
2014-5-1
We can't run your code, because you did not supply us with solarradiationcode. In which line do you expect there to be a branching?
My guess is that your problem is in one of your lines that is like
if N/3==i && i==2*N/3
You seem to be doing equality comparisons between floating-point numbers, which is not appropriate because not all numbers can be expressed exactly in floating-point arithmetic.
类别
在 帮助中心 和 File Exchange 中查找有关 PCM 的更多信息
2014-5-1
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
