Robin condition in a 1D FEM matlab code

Does anybody know any tutorial on implementing Robin conditions on a MATLAB 1D FEM code (which I did myself including creating NODES table, elemental stiffness matrix calculation and assembly etc). I am attaching the code for non-zero Dirichlet conditions in one boundary and free flow in the other (it is a advective-diffusive heat equation model). I am attaching also a link to the article whose results I am trying to replicate.

Thank you in advance. https://www.sciencedirect.com/science/article/pii/S0960148116303366 Code:

if true
  % code
%%Limpiamos todo
close all
clear all
clc
%%Introducimos los datos
NN = 2; % Vamos a trabajar con ELEMENTOS LINEALES (2 nodos por elemento)
ro = 1000; % kg/m3
Cp = 4185; % J/kg/K
k_pipe = 0.58; % W/m/K
v = 0.3; % m/s
r = 0.0131; % m
h = 1500; % W/m2/K
P = 2*pi*r; % m (Perímetro)
S = pi*r*r; % m^2 (Sección)
% La temperatura del terreno se extrae de otro documento o se toma cte
[Tg1,Tg2] = vector_Tg; % K
Long = 40;
hk = 1;
dt = 3;
t = 90;
alpha = k_pipe / (ro * Cp);
gamma = P * h / (ro * Cp * S);
Pe = v * hk / alpha; % Número de Petcler
beta = coth(Pe) - (1/Pe); %((exp(Pe)+exp(-Pe))/(exp(Pe)-exp(-Pe)))
% Definimos el número y posición de los nodos en x
COORDx = (Long:-hk:0)'; % Posiciones de los nodos de x
Mx = length(COORDx); % Número de nodos totales en x
NEx = Mx - 1; % Número de elementos
% Definimos el número y posición de los nodos en t
COORDt = (0:dt:max(t))';% Posiciones de los nodos de tiempo
Mt = length(COORDt); % Número de nodos totales de t o instantes de t
NEt = Mt -1; % Número de intervalos de tiempo
T = zeros(Mx,Mt); % Martiz que contendrá las soluciones
Tz = zeros(Mx,Mt); % Matriz que contendrá los valores de T desconocidos
Tw = zeros(Mx,Mt); % Matriz que contendrá los valores de T conocidos
% Componemos la tabla NODES
NODES = zeros(NEx,NN); ii = 0;
for i = 1:NEx
for j = 1:NN
ii = ii + 1;
NODES(i,j) = ii;
end
ii = ii - 1;
end
%%Condiciones de contorno
% Que en t=0, se sabe que T(z,0) = Tg(z,0)
T(:,1) = Tg1; % K
% Que la temperatura de entrada sea impuesta T.zin = Tin (Dirichlet)
Tw(1,:) = 300; % K
% Dado que el escape es de flujo libre (sin retorno, ni potencia añadida)
% la derivaba a la entrada y salida (los extremos) es 0.
B =0;
%%Cálculo de las matrices elementales
S = zeros(Mx); C = zeros(Mx); L1 = zeros(Mx);L2 = zeros(Mx);
[S_k,C_k,L_k] = Ak_Cte(hk,beta);
S_k = S_k * alpha;
C_k = C_k * v;
L1_k = L_k * gamma;
L2_k = L_k / dt;
%%Ensamblaje
for k=1:NEx
for i=1:NN
ig = NODES(k,i);
for j=1:NN
jg = NODES(k,j);
S(ig,jg) = S(ig,jg) + S_k(i,j);
C(ig,jg) = C(ig,jg) + C_k(i,j);
L1(ig,jg) = L1(ig,jg) + L1_k(i,j);
L2(ig,jg) = L2(ig,jg) + L2_k(i,j);
end
end
end
A = S + C + L1 + L2;
%%Resolvemos por el esquema de Euler implícito
% Para Tg constante
T(:,1) = Tg1;
for ti = 2:Mt
LD = L1*Tg1 + L2*T(:,ti-1) + B - A*Tw(:,ti);
Tz(2:Mx,ti) = A(2:Mx,2:Mx) \ LD(2:Mx);
T(:,ti) = [Tw(1,ti);Tz(2:Mx,ti)];
end
%%Graficamos la solución
figure(1), hold on, grid on
for i = 1:length(COORDt)
ylabel('z(m)'); xlabel('T(K)');
title({'Distribución de temperatura a lo largo de la tubería';'Tg = cte'});
plot(T(:,i),COORDx,'-');
axis([280 300 0 40]);
end
% Para Tg variable
T(:,1) = Tg2;
for ti = 2:Mt
LD = L1*Tg2 + L2*T(:,ti-1) + B - A*Tw(:,ti);
Tz(2:Mx,ti) = A(2:Mx,2:Mx) \ LD(2:Mx);
T(:,ti) = [Tw(1,ti);Tz(2:Mx,ti)];
end
%%Graficamos la solución
figure(2), hold on, grid on
for i = 1:length(COORDt)
ylabel('z(m)'); xlabel('T(K)');
title({'Distribución de temperatura a lo largo de la tubería';'Tg = Tg(z)'});
plot(T(:,i),COORDx,'-');
axis([280 300 0 40]);
function [S,C,L] = Ak_Cte(hk,beta)
%%Calculamos las componentes S, C y L de Ak
%%%Sk
S(1,1) = 1; S(1,2) = -1;
S(2,1) = -1; S(2,2) = 1;
S = 1/hk * S;
%%%Ck
C(1,1) = -1+beta; C(1,2) = 1-beta;
C(2,1) = -1-beta; C(2,2) = 1+beta;
C = 1/2 * C;
%%%Lk
L(1,1) = 4/6-3/6*beta; L(1,2) = 2/6-3/6*beta;
L(2,1) = 2/6+3/6*beta; L(2,2) = 4/6+3/6*beta;
L = hk/2 * L;
end
function [Tg1,Tg2] = vector_Tg
%%Si se toma Tg = constante
Tg1 = ones(41,1)*285;
%%Si se toma Tg = variable
Tg2=[10.4; 10.36; 10.32; 10.28; 10.24; 10.2; 10.16; 10.12; 10.08; 10.04;
10; 9.96; 9.92; 9.88; 9.84; 9.8; 9.76; 9.72; 9.68; 9.6; 9.61; 9.62;
9.64; 9.65; 9.67; 9.68; 9.69; 9.71; 9.72; 9.73; 9.75; 9.76; 9.77;
9.79; 9.8; 9.82; 8.9; 8; 7; 6; 5]+273.15*ones(41,1); Tg2 = fliplr(Tg2')';
end
end

Answer 1

[Neil Guertin]

You may find some useful information here:

https://www.mathworks.com/help/pde/ug/elliptic-pdes.html