finite volume method for 1D unsteady heat conduction

Hey All,

I am trying to simulate unsteady 1D heat conduction equation using MATLAB, I am following the instructions in the following link with changing one of the boundary conditions (West BC): https://www.youtube.com/watch?v=h3bHAEN7t6M. I also added heat to my equations, since I have heat generation for one of the three layers I have. The code runs normally for the first 500ish iterations, and then the plotted temperature suddenly disappears. Does this issue appear because of the values I'm feeding to the code or it is the convergence approach (lines 101-143)?

P.S

Even for the initial iterations, the temperature value appears insanely high

clear;close all;clc

Dvi=0.246926*2;            %ID of vessel (m)
Dvo=0.67*2;                %OD of vessel (m)
Dins=0.797*2;              %OD of insulation (m)
Dhyd=0.2;                  %hydraulic diameter of the coolant (air, in m)
L=1.4;                     %height of everything (m)

%Volume Descritizations [salt, vessel wall, insulator]
Discs=[50 30 10];
dx=[(0.5*Dvi)/Discs(1)...
    (0.5*(Dvo-Dvi))/Discs(2)...
    (0.5*(Dins-Dvo))/Discs(3)];        %Spatial increment (dx)

Ac= 0.25.*pi*dx.^2;                    %Cross-sectinal area

%Time descritization
dt=10;                                 %time increment (sec)
t=1:dt:1000;                           %time vector (sec)

%Heat rate
q=120975.*t.^-0.293;

%Thermal properties of salt
cps=4683718.0;                         %specific heat (J/m3-K)
rhov=8e3;                              %density (kg/m3)
cpv=490;                               %heat capacity for SS304
cpins=144640;                          %specific heat (J/m3-K)

%Thermal conductivities
kins=0.15;                             %W/m-K
kv=14;                                 %W/m-K
ks=20;                                 %W/m-K

%Coolant (air) temperature initialization
T0 =100;                               %T initial
Tinf=300;                              %temperature of air (K)
hinf=12;                               %assumed (W/m2-K)
Rinf=1/(hinf*pi*Dhyd*L);               %Thermal resistance due to convection in coolant
K_a     = 26.24;                       %Thermal conductivity of the air
Cp_a    = 1e3;                         %heat capacity for the air
rho_a   = 1.184;                       %density of the air
mu_a    = 18.6e-6;                     %dynamic viscosity of the air
nu      = mu_a/rho_a;
pr_a    = mu_a*Cp_a/K_a;               %prandtl number for the salt
g       = 9.81;
alp=K_a/rho_a*Cp_a;
beta=3.32e-3;                          %expansion coefficient puuled out of the reference


%% Constructing the A matrix
A=eye(sum(Discs));
%%%%%%%%%%%%%%%%INPUTS

Adiag_ins=-(2*kins*Ac(3)./dx(3) + cpins*Ac(3)*dx(3)./dt);
A_offDiag_ins=kins*Ac(3)./dx(3);


Adiag_v=-(2*kv*Ac(2)./dx(2) + rhov*cpv*Ac(2)*dx(2)./dt);
A_offDiag_v=kv*Ac(2)./dx(2);


Adiag_s=-(2*ks*Ac(1)./dx(1) + cps*Ac(1)*dx(1)./dt);
A_offDiag_s=ks*Ac(1)./dx(1);
%%%%%%%%%%%%%%%%

for i=1:Discs(3)
    A(i,i)=Adiag_ins;
    [A(i,i+1),A(i+1,i)]=deal(A_offDiag_ins);
end

% A(Discs(3),Discs(3)-1)=A_offDiag_v;

for i=Discs(3):sum([Discs(3) Discs(2)])
    A(i,i)=Adiag_v;
    [A(i,i+1),A(i+1,i)]=deal(A_offDiag_v);
end

% A(sum([Discs(3) Discs(2)]),sum([Discs(3) Discs(2)])-1)=A_offDiag_s;

for i=sum([Discs(3) Discs(2)]):sum(Discs)-1
    A(i,i)=Adiag_s;
    [A(i,i+1),A(i+1,i)]=deal(A_offDiag_s);
end

% A(i+1,i+1)=Adiag_s;
A=rot90(A,2);
A(1,1)=-(ks*Ac(1)./dx(1) + cps*Ac(1)*dx(1)./(2*dt));
A(end,end)=-(kins*Ac(3)./dx(3) + cpins*Ac(3)*dx(3)./(2*dt));

% imagesc(A); colormap jet

%% Solving for temperature
T=zeros(sum(Discs),length(t))+Tinf;

ResidualTolerance=0.1;

close all;clc
figure('Units','inches','Position',[1 1 10 5]); hold on
for i=1:length(t)-1
    ct=0; TMax = 500;
    while abs(max(T(:,i+1))-TMax)>ResidualTolerance
        ct=ct+1;
        TMax=max(T(:,i+1));
        
        %%%%%%%%%%%%%%%%INPUTS
        B_s=-q(i)-(cps*Ac(1)*dx(1)./dt);
        B_v=-(rhov*cpv*Ac(2)*dx(2)./dt);
        B_ins=-(cpins*Ac(3)*dx(3)./dt);
        %%%%%%%%%%%%%%%%BMatrix
        BConst=[repmat(B_s,[Discs(1),1]) ; ...
            repmat(B_v,[Discs(2),1]) ; ...
            repmat(B_ins,[Discs(3),1])];
        
        
        %%%%%%%%%%Boundaries for vector b
        B_boundary1=-q(i)-((cps*Ac(1)*dx(1)./(2*dt))*TMax);
        B_boundary2=-((cpins*Ac(3)*dx(3)./(2*dt))*T(end,i))-(hinf*Ac(3)*Tinf);
        
        
        B(:,i)=T(:,i).*BConst;
        B(1,i)=B_boundary1;
        B(end,i)=B_boundary2;
%         C = inv(A);
        T(:,i+1)=A\B(:,i);
        
        Ra=(g.*beta./(nu*alp)).*(Tinf-T(end,i+1)).*(Dhyd^3); %initial guess for Ra
        Nu=(0.825+(0.387*Ra.^(1/6))/(1+(0.492/pr_a).^(9/16)).^(8/27)).^2;
        h=Nu.*K_a./Dhyd;
        Rinf=1/(hinf*pi*Dhyd*L); %Thermal resistance due to convection in coolant
        Tinf=-(Rinf*q(i))+T(end,i+1);
        
    end
    
    disp([i ct])
    cla;
    plot(T(:,i+1),'k');
    ylabel('$T$','Interpreter','latex');
    pause(0.1)
    title(num2str(t(i+1)));
    
end
hold off;

[Sargondjani]

You can not expect us to check every line of your code. It is not even clear whether you suspect the formulas are wrong, or if there is some coding issue.

You should be more precise where the problem might occur: find out yourself where the problem exactly is (or might be).

[Laith Zaidan]

I think if you read the description carefully, you will see that I have already mentioned the lines I'm suspicious about. My problem is more of coding logic, not coding errors.

[KSSV]

First check would are all the constants or variables in the same unit?

Answer 1

[Walter Roberson]

On iteration 52,

        T(:,i+1)=A\B(:,i);

generates infinities. It involves quantities that would become about 2.55e+308 which is larger than the largest representable number of about 1.8E+308 .

A is full rank and a reasonable condition number at that point, but B at that point involves values as small as -1.1E+307

Once you have infinity generated, then most of the time either the next array operation or the one after that will generate NaN due to attempts to take the difference of infinities.

Plots do not draw nan or infinite values.

[Laith Zaidan]

Alright, thank you very much