Implicit heat conduction 1D

Question

0 个投票

I am trying to simulate 1d heat transfer thru wall of 3mm thickness. One one side I have convection,radiation and solar heat flux. On other side I have fluid at 20K which is in contact with wall. I have discretized whole domain and written matrices and tried solving for temperatures. I am getting some values but that are not correct.

i have attached the derivation also.

dt = 0.01;                      % time step size
time = 150;                     % time before launch
nt = time/dt + 1 ;                  % number of time steps
t = linspace(0,time,nt);        % time array
L = 3e-3;
nx = 10;
x = linspace(0,L,nx);
dx = x(2)-x(1);
rho_w = 7850;
k_l = 0.101 ;        % LH2 thermal conductivity in W/m-K
C_l = 10.019e3;         % Specific heat capacity (J/kg-K) 
rho_l = 70.35;
Cw_range = [13.65, 275,360,425]; % Specific heat capacity (J/kg-K) corresponding to temperature range
T_Cp = [20,100,150,200];        % Temperatures (K) for specific heat capacity
h_air = 4.3;         % Heat transfer coefficient for air (W/m^2-K)
T_amb = 300;         % Ambient air temperature (K)
e_c   = 0.3;         % Emissivity of Inconel
sigma = 5.67e-8;     % Stefan-Boltzmann constant (W/m^2-K^4)
sol = 1400;
T_LH2 = 20;          % LH2 temperature (K)
% Variable thermal conductivity
p1 = -3.0277e-9;   
p2 = 2.7759e-6;   
p3 = -9.4393e-4;  
p4 = 0.1682;   
p5 = -0.6885;
% k_w = p1 * T_old_steady.^4 + p2 * T_old_steady.^3 + p3 * T_old_steady.^2 + p4 * T_old_steady + p5;
T_range = linspace(20,1000,500);
k_values = p1 * T_range.^4 + p2 * T_range.^3 + p3 * T_range.^2 + p4 * T_range + p5;
k_values(T_range > 300) = p1 * 300^4 + p2 * 300^3 + p3 * 300^2 + p4 * 300 + p5;
% Precompute specific heat capacity (Cw) over the temperature range
Cw_values = interp1(T_Cp, Cw_range, T_range, 'linear', 'extrap');
Cw_values(T_range >= 300)=interp1(T_Cp, Cw_range, 300, 'linear', 'extrap');
% Initialize temperature
T_initial = ones(nx+1,1) * T_amb;
T_initial(nx+1) = T_LH2;
T_old = T_initial;
% Create a matrix to store temperature at each time and location
% T_all = zeros(nt, nx+1);  % nt rows, nx+1 columns
for n = 1 : nt
    k_w = interp1(T_range, k_values, T_old, 'linear', 'extrap');
    C_w = interp1(T_range, Cw_values, T_old, 'linear', 'extrap');
    
    alpha_w = k_w ./ (C_w * rho_w);
    
    Alpha = (k_w*dt) ./ (rho_w*C_w*dx^2);
    alpha_f = (k_l*dt)/(rho_l*C_l*dx^2);
       
    A = zeros(nx+1,nx+1);
    B = zeros(nx+1,1);
    for i = 2:nx-1
    alpha_left = (2*Alpha(i)*Alpha(i-1))/(Alpha(i) + Alpha(i-1));
    alpha_right = (2*Alpha(i)*Alpha(i+1))/(Alpha(i) + Alpha(i+1));
 
    
    A(i,i-1) = -alpha_left;
    A(i,i) = 1+alpha_right+alpha_left;
    A(i,i+1) = -alpha_right;
    end
      
      % Left boundary
      
      k_eff_left = interp1(T_range, k_values, T_initial(1), 'linear', 'extrap');
      Cw_left = interp1(T_range, Cw_values, T_initial(1), 'linear', 'extrap');
      coeff = (rho_w * Cw_left * dx) ./ (2 * dt);
      K1 = coeff*T_old(1) - (h_air*T_amb) + (3*e_c*sigma*(T_old(1))^4) + (e_c*sigma*T_amb^4);
      A(1,1) = coeff - h_air + 4*e_c*sigma*(T_old(1))^3 + (2*k_eff_left/dx);
      A(1,2) = -2*k_eff_left/dx;
      B(1) = K1;
      A(nx,nx) =  1 + 8 * alpha_right;
      A(nx,nx-1) = - 4 * alpha_right;
      A(nx,nx+1) = - 4 * alpha_right;
      B(nx) = T_old(nx);
      % Right boundary (connection to LH2 tank)
      k_eff_right = interp1(T_range, k_values, T_initial(end), 'linear', 'extrap');
      A(nx+1,nx) = -alpha_f;
      A(nx+1,nx+1) = 1 + alpha_f;
      B(nx+1) = T_LH2;
      
      a = diag(A);      % main diagonal 
      b = diag(A,1);    % upper diagonal 
      c = diag(A,-1);   % lower diagonal 
    
    % Forward substitution
      for i = 2:nx
          factor = c(i-1) / a(i-1);
          a(i) = a(i) - factor * b(i-1);
          B(i) = B(i) - factor * B(i-1);
      end
  % Backward substitution
  T = zeros(nx+1,1);
  T(nx+1) = B(nx+1) / a(nx+1);
    for i = nx:-1:1
        T(i) = (B(i) - b(i) * T(i+1)) / a(i);
    end
    T_initial = T;
    T_old = T_initial; 
%     T_all(n, :) = T';
end