help with implicit method

1 次查看（过去 30 天）
Libya
Libya 2014-5-23
编辑： Libya 2014-5-25
Hi guys, This is my first code that I have ever written. I have wrote this code, and I took almost three months to get to this stage. the code is running good but the results are very high, where the output temperature of the material is 250C. this code has a subroutine code where it is consider as an input (Gt). I have changed all the material properties, but nothing change in the output, when I divided Gt (subroutine code)by 10, the output temperature has been divided by 10 and became 25 C, which is acceptable. there is nothing to do with the subroutine, because it is absolutely right.
Is there anybody could help with that please, and find out why my output temperature is very high.
clear all
clc
ho=50;
hi=10;
N=27; %Number of space intervals
M=3600; %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=2.435;
dT=0.2;
Tm=5.5;
Gt=solarradiationcode;
dz=2;
xr = length(Gt);
Te = zeros(xr,1);
Final(1,xr) = zeros;
dt =3;
%initial condition
u(:,1)=ones(N,1)*15;
for h = 1:length(Gt)
Te(h) = To*(abs*Gt(h)/ho);
for j=1:1:M %time step
alfa= k_b/(Cp_b*roh_b);
% define the ratio r
r=alfa*dt/dx^2; %Fourier number
Beta=[2*r*bo*Te(h); 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:1:N-1
%
%
if i<N/3 || i>2*N/3
'1'
alfa= k_b/(Cp_b*roh_b);
r=alfa*dt/dx^2;
end
%
if i==N/3 || i==2*N/3
% % '2'
%
Cp=2*Cp_b*Cp_ps/(Cp_ps+Cp_b);
roh=2*roh_ps*roh_b/(roh_ps+roh_b);
k=2*k_b*k_ps/(k_ps+k_b);
alfa= k/(Cp*roh);
r=alfa*dt/dx^2;
end
if i>N/3 && i<2*N/3
% '3'
T=u(i,j);
if T<Tm
Cp=Cp_ps;
roh=roh_ps;
k=k_ps;
elseif T>Tm-dT && T<Tm+dT
Cp= ((Cp_ps+Cp_pl)/2)+((roh_ps+roh_pl)/2)*(lamda/dT);
roh=roh_ps;
k=k_ps;
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,UU] = LU(A);
% note: matlab can tell the difference between u and U, so I could have U as the
% matrix and u as the solution. But not all languages will do this, sometimes you
% have to tell them explicitly to care about case. The safest thing is to not have the
% habit of depending on cases.
% then we'll just do
%[y] = down_solve(L,u_old);
%[u_new] = up_solve(UU,y);
b=u(:,j)+Beta;
[y] = down_solve(L,b);
u(:,j+1) = up_solve(UU,y);
end
Final(1,h) = u(15,M);
% Final(1,h) = u(15,M);
%Final(1,h) = u(15,M);
end
plot(Final)
  4 个评论
显示 2更早的评论
Image Analyst
Image Analyst 2014-5-24
Well, attach them. Attach all needed m-files with the paper clip icon so we can try your code. Also, read this: http://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup
Libya
Libya 2014-5-25
编辑：Libya 2014-5-25
When I went through this code line by line I found that, when I eliminate the A matrices at each time step the output increases. Anybody could help.

请先登录，再进行评论。

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

回答（0 个）

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by