Interface boundary condition in transient heat problem

Hello,

It is my first problem and first post on this forum. I am trying to solve 1D transient heat equation in multilayer electrical cable using Crank-Nicolson method, and I have problems with interface boundary conditions. I'm not sure how to apply them. I'm at very begining of learning so my code is far away from ideal, but at this point it is important for me to resolve my problem. I will optimize the code later. Here is my code:

clc; clear;

% -------- Simulation Parameters -------- %

L_i = 0.021;                            % Insulation + inner and outer screen thickness
L_ls = 0.0035;                          % lead sheat thickness
L_PE = 0.0035;                          % PE sheat thickness
L_st = 0.0005;                          % steal tape thickness
L_b = 0.0005;                           % bedding thickness
L_ta = 0.0055;                          % armor thickness
L_j = 0.004;                            % jute thickness
L_g = 0.0385;                           % ground thickness

N  = 154;
Lx = 0.077;
dx = Lx/N;
x  = 0:dx:Lx;

M = 1420;
tf = 1420;     
dt = tf/M;
t = 0:dt:tf;

% -------- thermal diffusivity -------- %

alpha_c = 1.15e-4;                      % conductor
alpha_cs = 9.83e-8;                     % conductor screen
alpha_i = 9.83e-8;                      % insulation
alpha_is = 9.83e-8;                     % insulation screen
alpha_ls = 2.32e-5;                     % lead sheat
alpha_PE = 1.18e-7;                     % PE sheath
alpha_st = 1.64e-5;                     % steel tape
alpha_b = 3.18e-7;                      % bedding
alpha_ta = 1.64e-5;                     % tensile armour
alpha_j = 8.79e-8;                      % jute
alpha_g = 5e-7;                         % ground

% -------- diffusion number -------- %

d_c = (alpha_c*dt)/(dx^2);              % conductor
d_cs = (alpha_cs*dt)/(dx^2);            % conductor screen
d_i = (alpha_i*dt)/(dx^2);              % insulation
d_is = (alpha_is*dt)/(dx^2);            % insulation screen
d_ls = (alpha_ls*dt)/(dx^2);            % lead sheath
d_PE = (alpha_PE*dt)/(dx^2);            % PE sheath
d_st = (alpha_st*dt)/(dx^2);            % steel tape
d_b = (alpha_b*dt)/(dx^2);              % bedding
d_ta = (alpha_ta*dt)/(dx^2);            % tensile armour
d_j = (alpha_j*dt)/(dx^2);              % jute
d_g = (alpha_g*dt)/(dx^2);              % ground

% -------- thermal conductivity -------- %
k_i = 0.2;
k_ls = 33;
k_PE = 0.3;
k_st = 60;
k_b = 0.15;
k_ta = 59;
k_j = 0.18;

% -------- other input data -------- %
r_c = 0.0265;                           % conductor radius
I = 1400;                               % current in conductor
R = 8.19e-06;                           % conductor resistance
q = (I^2 * R)/(2 * pi * r_c);           % heat flux
Tin = 300;                              % Starting temperature
T_out = Tin;                            % outside temperature

% -------- Matrix definition -------- %

A = zeros(N,N);
B = zeros(N,N);
dT = zeros(N,M);
BC = zeros(N,1);                        % Boundary conditions
IC = @(x) Tin;                          % Initial Conditions
dT(:,1) = IC(x);                        % apply of initial condition


% ------------------- A Matrix ------------------- %

for j = 1:N+1                               % Column loop
    for i = 1:(L_i/dx)                      % Row loop for insulation
        if i==j && i>1
            A(i,j) = 2*(1+d_i);
        elseif i == 1 && i == j
            A(i,j) = 2*(1+d_i);
        elseif i + 1==j && j>2
            A(i,j) = -d_i;
        elseif i + 1==j && j<=2
            A(i,j) = -2*d_i;
        elseif i - 1 == j
            A(i,j) = -d_i;
        end
    end
    for i = (L_i/dx)+1:(L_i+L_ls)/dx        % Row loop for lead sheat
        if i==j
             A(i,j) = 2*(1+d_ls);
        elseif i==j+1
            A(i,j) = -d_ls;
        elseif i==j-1
            A(i,j) = -d_ls;
         end
    end
    for i = (L_i + L_ls)/dx +1 : (L_i + L_ls + L_PE)/dx        % Row loop for PE sheat
        if i==j
            A(i,j) = 2*(1+d_PE);
        elseif i==j+1
            A(i,j) = -d_PE;
        elseif i==j-1
            A(i,j) = -d_PE;
        else
         end
    end
    for i = (L_i + L_ls + L_PE)/dx +1 : (L_i + L_ls + L_PE + L_st)/dx     % Row loop for steel tape
        if i==j
            A(i,j) = 2*(1+d_st);
        elseif i==j+1
            A(i,j) = -d_st;
        elseif i==j-1
            A(i,j) = -d_st;
        else
         end
    end
    for i = (L_i + L_ls + L_PE + L_st)/dx +1 : (L_i + L_ls + L_PE + L_st + L_b)/dx     % Row loop for bedding
        if i==j
            A(i,j) = 2*(1+d_b);
        elseif i==j+1
            A(i,j) = -d_b;
        elseif i==j-1
            A(i,j) = -d_b;
         end
    end
    for i = (L_i + L_ls + L_PE + L_st + L_b)/dx +1: (L_i + L_ls + L_PE + L_st + L_b + L_ta)/dx     % Row loop for armour
        if i==j
            A(i,j) = 2*(1+d_ta);
        elseif i==j+1
            A(i,j) = -d_ta;
        elseif i==j-1
            A(i,j) = -d_ta;
         end
    end
    for i = (L_i + L_ls + L_PE + L_st + L_b + L_ta)/dx +1 : (L_i + L_ls + L_PE + L_st + L_b + L_ta + L_j)/dx       % Row loop for propylen and bitumen
        if i==j
            A(i,j) = 2*(1+d_j);
        elseif i==j+1
            A(i,j) = -d_j;
        elseif i==j-1
            A(i,j) = -d_j;
        end
    end
    for i = 78 : N+1                    % Row loop for ground
        if i == j && i <= N-1
           A(i,j) = 2*(1+d_g);
        elseif i == j+1 && i <= N-1
            A(i,j) = -d_g;
        elseif i == j-1 && i <= N-1
            A(i,j) = -d_g;
        else
            A(154,154) = 1;
        end
    end
end

% ------------------- B Matrix ------------------- %
for j = 1:N+1
    for i = 1:(L_i/dx)
        if i==j && i>1
            B(i,j) = 2*(1-d_i);
        elseif i==j+1
            B(i,j) = d_i;
        elseif i + 1==j && j<=2
            B(i,j) = 2*d_i;
        elseif i==j-1 && j>2
            B(i,j) = d_i;
        elseif i==1 && j==1
            B(i,j) = 2*(1-d_i); %%1;
        end
    end
    for i = (L_i/dx)+2:(L_i+L_ls)/dx +1
        if i==j && i>1
            B(i,j) = 2*(1-d_ls);
        elseif i==j+1
            B(i,j) = d_ls;
        elseif i==j-1
            B(i,j) = d_ls;
         end
    end
    for i = (L_i + L_ls)/dx +1 : (L_i + L_ls + L_PE)/dx 
        if i==j
            B(i,j) = 2*(1-d_PE);
        elseif i==j+1
            B(i,j) = d_PE;
        elseif i==j-1
            B(i,j) = d_PE;
         end
    end
    for i = (L_i + L_ls + L_PE)/dx +1 : (L_i + L_ls + L_PE + L_st)/dx 
        if i==j
            B(i,j) = 2*(1-d_st);
        elseif i==j+1
            B(i,j) = d_st;
        elseif i==j-1
            B(i,j) = d_st;
         end
    end
    for i = (L_i + L_ls + L_PE + L_st)/dx +1 : (L_i + L_ls + L_PE + L_st + L_b)/dx 
        if i==j
            B(i,j) = 2*(1-d_b);
        elseif i==j+1
            B(i,j) = d_b;
        elseif i==j-1
            B(i,j) = d_b;
         end
    end
    for i = (L_i + L_ls + L_PE + L_st + L_b)/dx +1: (L_i + L_ls + L_PE + L_st + L_b + L_ta)/dx 
        if i==j
            B(i,j) = 2*(1-d_ta);
        elseif i==j+1
            B(i,j) = d_ta;
        elseif i==j-1
            B(i,j) = d_ta;
         end
    end
    for i = (L_i + L_ls + L_PE + L_st + L_b + L_ta)/dx +1 : (L_i + L_ls + L_PE + L_st + L_b + L_ta + L_j)/dx 
        if i==j
            B(i,j) = 2*(1-d_j);
        elseif i==j+1
            B(i,j) = d_j;
        elseif i==j-1
            B(i,j) = d_j;
         end
    end
    for i = 78 : N+1
        if i==j && i<= N-1
            B(i,j) = 2*(1-d_g);
        elseif i==j+1 && i <= N-1
            B(i,j) = d_g;
        elseif i==j-1 && i <= N-1
            B(i,j) = d_g;
        else
            B(154,154) = 1;
        end
    end
end

for j = 1:N                     % Boundary condition loop
    if j ==1
       BC(j,1) = 4 * d_i * (q/k_i) * dx;
    else
        BC(j,1) = 0;
    end
end

for j = 2:M               % Time loop
    for i = 1:N           % Space loop
        dT(:,j) = A\(B*dT(:,j-1)+BC);
    end
end

I hope this code is ok at this stage.

In row T0 there is boundary from constant heat flux from conductor (I use a centered difference approximation and put it into T0) and in node 154 is a fixed soil temperature. I will be very grateful if someone could show me (in a clear way :) ) how to add this boundaries. I'm not sure but I think it should be something from centered difference approximation like T42(-k1-k2) = -k2T44 - k1T40?

I work on this code in Matlab without any toolboxes.

Thank you all in advance for helping me with above.

Best regards

Kamil

[Torsten]

Use "pdepe" and choose the interface coordinates as grid points.

This will work if the transient conditions are continuity of temperature and heat flux.

If you want to put jump discontinuities, the problem will become more difficult.

[Kamil]

Thank you for your reply Torsten, I really appreciate it. My goal is to apply a steady stream of heat from the current for eight hours, then remove the current and cool the cable for a net sixteen hours. I want to compare my model with the factory tests of the cable.

[Torsten]

And why do you need to consider interfaces (i.e. intermediate points between start and end point of the cable) for this problem ?

[Kamil]

Because I mainly want to plot the distribution of the electric field in the insulation. Since this is a DC cable, the electric field depends on the temperature distribution. I would like to study how it behaves under different operating conditions of the cable.

[Torsten]

I'm too lazy to extract the interface conditions from your code above.

As said, they usually are

T(xi-) =  T(xi+)
lambda_l * dT/dx (xi-) = lambda_r * dT/dx (xi+)

which means: continuity of temperature and heat flux at the interface.

If you give your two conditions at the interfaces in a similar manner, I'll try to tell you about their discretization.

[Kamil]

The continuity of temperature and heat flux at the interface of the layers is exactly what I wanted to do. Since I have the same grid space for the whole cable, I thought it should be like that: T(i+1) = T(i-1). When I tried to apply this boundary condition, I got graphs like the first and second (first is a temperature at first point of insulation and a second one is a temperature along of cable after 24 hours). Without the boundary condition, I got graphs like the third and fourth. I doing something wrong in applying those boundary conditions and that is why I asking for help. Graphs 3 and 4 are close to the truth, but the temperature drop along the cable is too high.

[Torsten]

If your equations are

rho_l*cp_l * dT/dt = lambda_l * d^2T/dx^2, x <= xi

with discretization points

x_l1 < x_l2 < ... < x_ln = xi

and

rho_r*cp_r * dT/dt = lambda_r * d^2T/dx^2,  x >= xi

with discretization points

x_r1 = xi < x_r2 < x_r3 < ... < x_rm

the usual way to incorporate the interface condition is to introduce ghost points

x_l_ghost = xi + (xi - x_l(n-1))
x_r_ghost = xi - (x_r2 - xi)

discretize the interface condition as

lambda_l * dT/dx (xi-) = lambda_l * (T_l_ghost - T_l(n-1)))/(x_l_ghost - x_l(n-1))     (1)
lambda_r * dT/dx (xi+) = lambda_r * (T_r2 - T_r_ghost)/(x_r2 - x_r_ghost)              (2)

and discretize the two PDEs in x = xi as

rho_l*cp_l * dT/dt (xi-) = lambda_l*((T_l_ghost - Ti)/(x_l_ghost - xi) - (Ti - T_l(n-1))/(xi - x_l(n-1)))/((x_l_ghost+xi)/2 - (xi+x_l(n-1))/2)      (3)
rho_r*cp_r * dT/dt (xi+) = lambda_r*((T_r2) - Ti)/(x_r2 - xi) - (Ti - T_r_ghost)/(xi - x_r_ghost))/((x_r2+xi)/2 - (xi+x_r_ghost)/2)                 (4)

Since the temperatures at the interfaces agree, you arrive at the following three equations for Ti, T_l_ghost and T_r_ghost:

(a) : right-hand side of equation (1) == right-hand side of equation (2) (continuity of heat flux)

(b) : right-hand-side of equation (3) / (rho_l*cp_l) == right-hand-side of equation (4) / (rho_r*cp_r) (continuity of temperature)

(c) : one of the differential equations (3) or (4) to compute the temperature at the interface.

The easiest way to solve this system is to solve the linear algebraic equations (a) and (b) for T_l_ghost and T_r_ghost and insert the obtained expressions in the differential equation (c). (c) will then be the differential equation for Ti and will only depend on known quantities.

Answer 1

[Torsten]

Here is a sample code for what I wrote above:

  % Set parameters
  x1start = 0.0;
  x1end = 0.25;
  x2start = x1end;
  x2end = 1.0;
  nx1 = 50;
  nx2 = 150;
  x1 = linspace(x1start,x1end,nx1).';
  x2 = linspace(x2start,x2end,nx2).';
  rho1 = 100;
  cp1 = 10;
  lambda1 = 0.4;
  rho2 = 20;
  cp2 = 25;
  lambda2 = 10;
  
  % Set initial conditions
  Tstart = 200;
  Tend = 400;
  T0 = 0;

  nodes = nx1+nx2-1;
  y0 = zeros(nodes,1);
  y0(1) = Tstart;
  y0(2:nx1+nx2-2) = T0;
  y0(nx1+nx2-1) = Tend;

  % Set time span of integration
  tspan = 0:10:100;

  [T,Y] = ode15s(@(t,y)fun(t,y,x1,nx1,x2,nx2,lambda1,rho1,cp1,lambda2,rho2,cp2),tspan,y0);
  
  % Plot results 
  x = [x1;x2(2:end)];
  plot(x,Y.')

  function dy = fun(t,y,x1,nx1,x2,nx2,lambda1,rho1,cp1,lambda2,rho2,cp2)
    T1 = y(1:nx1);          % Temperature zone 1
    T2 = y(nx1:nx1+nx2-1);  % Temperature zone 2

    % Compute temperature in ghost points    
    xg1 = x1(end)+(x1(end)-x1(end-1));
    xg2 = x2(1)-(x2(2)-x2(1));
    % Set up linear system of equations for [Tg1, Tg2]
    a11 = lambda1/(xg1-x1(nx1-1));
    a12 = lambda2/(x2(2)-xg2);
    b1 = lambda1*T1(nx1-1)/(xg1-x1(nx1-1))+lambda2*T2(2)/(x2(2)-xg2);

    a21 = (lambda1/(xg1-x1(nx1))/((x1(nx1)+xg1)/2 - (x1(nx1)+x1(nx1-1))/2))*rho2*cp2;
    a22 = (-lambda2/(x2(1)-xg2)/((x2(2)+x2(1))/2 - (x2(1)+xg2)/2))*rho1*cp1;
    b2 = (lambda1*T1(nx1)/(xg1-x1(nx1))+lambda1*(T1(nx1)-T1(nx1-1))/(x1(nx1)-x1(nx1-1)))/...
          ((x1(nx1)+xg1)/2-(x1(nx1)+x1(nx1-1))/2)*rho2*cp2 +...
         (lambda2*(T2(2)-T2(1))/(x2(2)-x2(1))-lambda2*T2(1)/(x2(1)-xg2))/...
          ((x2(2)+x2(1))/2 - (x2(1)+xg2)/2)*rho1*cp1;

    A = [a11,a12;a21,a22];
    b = [b1;b2];
    % Solve linear system for [Tg1, Tg2]
    sol = A\b;

    Tg1 = sol(1);
    Tg2 = sol(2);

    % Solve heat equation in the two zones
    % Zone 1
    T1 = [T1;Tg1];
    x1 = [x1;xg1];
    dT1dt(1) = 0.0;
    for i=2:numel(x1)-1
      dT1dt(i) = (lambda1*(T1(i+1)-T1(i))/(x1(i+1)-x1(i)) -...
                  lambda1*(T1(i)-T1(i-1))/(x1(i)-x1(i-1)))/...
                  ((x1(i+1)+x1(i))/2-(x1(i)+x1(i-1))/2)/(rho1*cp1);
    end
    % Zone 2
    for i=2:numel(x2)-1
      dT2dt(i) = (lambda2*(T2(i+1)-T2(i))/(x2(i+1)-x2(i)) -...
                  lambda2*(T2(i)-T2(i-1))/(x2(i)-x2(i-1)))/...
                   ((x2(i+1)+x2(i))/2-(x2(i)+x2(i-1))/2)/(rho2*cp2);
    end
    dT2dt(end+1) = 0.0;
    
    % Return time derivatives
    dy = [dT1dt(1:end),dT2dt(2:end)];
    dy = dy.';
  end

[Kamil]

I really appreciated your effort. I believe it will help me reorganize my code to make it work properly. Thanks a lot Torsten!

Best Regards

Kamil