2D steady state thermal model using Volume Element Method (VEM) for a cable with flowing fluid in center

Hello,

I am modelling the temperture distriubtion in a superconduting cable via the volume element method. See the image attached for a describtion of the VEM thermal network.

The temperture values are calcualted at every layer of the cable using the tempeerture values of all the neighbouring elements. The cable is descriitzed into a matrix of 7xn elements where n is the max number of elements along the length of the cable. below is my code, however the reuslts i recieve do not make sense as one would expect a gradual increase in temperure along the leneght of the cable, howver the results show that the temperure jumps rapidly and decreases depending on the layer.

See here

clc;
clear all;

% Cable Dimensions
L = 200000; % Length of the cable
N = 100; % Number of nodes
dz = L/N; % Length of each element

z = 0:dz:L;

% Initialize temperature values
T1 = ones(N+1, 1)*14; % Temperature values at each node

T2 = ones(N+1, 1)*14;

T3 = ones(N+1, 1)*14;

T4 = ones(N+1, 1)*14;

T5 = ones(N+1, 1)*20;

T7 = ones(N+1, 1)*277;

T8 = ones(N+1, 1)*278;


T = [T1,T2,T3,T4,T5,T7,T8];

T(1,1)=14;
T(1,2)=14;
T(1,3)=14;
T(1,4)=14;
T(1,5)=20;
T(1,6)=277;
T(1,7)=278;



%Simuntaneosuly solving array of 7XN equations

options = optimoptions(@fsolve,'MaxFunctionEvaluations', 5000, 'MaxIterations', 5000);

[T_solved,fval] = fsolve(@energy_balance,T,options);

%Plotting Temperatures[1-7] as afunction of z

figure (1);
plot(T_solved(:,1));
title('T1 as  afunction of cable distance z');
xlabel('number of nodes'), ylabel('T1 (K)');
figure (2);
plot(T_solved(:,2));
title('T2 as  afunction of cable distance z');
xlabel('number of nodes'), ylabel('T2 (K)');
figure (3);
plot(T_solved(:,3));
title('T3 as  afunction of cable distance z');
xlabel('number of nodes'), ylabel('T3 (K)');
figure (4);
plot(T_solved(:,4));
title('T4 as  afunction of cable distance z');
xlabel('number of nodes'), ylabel('T4 (K)');
figure (5);
plot(T_solved(:,5));
title('T5 as  afunction of cable distance z');
xlabel('number of nodes'), ylabel('T5 (K)');
figure (6);
plot(T_solved(:,6));
title('T7 as  afunction of cable distance z');
xlabel('number of nodes'), ylabel('T7 (K)');
figure (7);
plot(T_solved(:,7));
title('T8 as  afunction of cable distance z');
xlabel('number of nodes'), ylabel('T8 (K)');

%Plotting Residauls[F1-F7] 

figure (8);
hold on;
for c=1:size(fval,2)
    plot(fval(:,c));
end
hold off;
title('Residuals at each element');
xlabel('number of nodes'), ylabel('Residuals');
legend('F1', 'F2', 'F3', 'F4', 'F5', 'F6', 'F7');

% ---------------------Equation Balance function---------------------------------------------
function F = energy_balance(T)


% Cable Dimensions
L = 200000; % Length of the cable
N = 100; % Number of nodes
dz = L/N; % Length of each element

T_amb0 = 278; %Ambient Sea Temp
A_pipeInt=0.15707963;
% Heat Transfer constants
h_conv = 0.78;
sigma = 5.67E-08;
ep_5 = 0.05;
ep_7 = 0.05;
cp_int = 10;
mdot_int = 0.01663667;


% Energy Equations Constants
E2A = 3.1895E-06;
E2B = 0.000596023;

E3A = 0.016648776;
E3B = 3.18951E-06;
E3C = 1.0032E-07;
Q_gen= 0.001;

E4A = 0.00051836;
E4B = 1.0032E-07;
E4C = 0.000148357;

E5A = 3.40234E-07;
E5B = 0.000148357;

A56 = 18.22124;
A67 = 81.0553;
F22 = 0.2248;
F21 = 0.7752;

E7A = 0.148911492;
E7B = 1.1955E-05;

E8A = 2.0656E-05;
E8B = 1.1955E-05;
E8C = 1.18111E-05;

%intialie equation matrix F
F =zeros(7*N,1);

    for i = 1:N

       if i==1
           %Initial tempertures
T(1,1)=14;
T(1,2)=14;
T(1,3)=14;
T(1,4)=14;
T(1,5)=20;
T(1,6)=277;
T(1,7)=278;

       else
    

    
    F(i,1) = h_conv*A_pipeInt*(T(i,2)-T(i,1)) + mdot_int*cp_int*(T(i-1,1)-T(i,1));

    F(i,2) = -h_conv*A_pipeInt*(T(i,2)-T(i,1)) + E2A*(T(i,3)-T(i,2)) + E2B*(T(i,2)-T(i-1,2)) + E2B*(T(i+1,2)-T(i,2));

    F(i,3) = -E3A*(T(i,3)-T(i-1,3)) + E3A*(T(i+1,3)-T(i,3)) - E3B*(T(i,3)-T(i,2)) - E3C*(T(i,3) -T(i,4)) + Q_gen;

    F(i,4) = -E4A*(T(i,4)-T(i-1,4)) + E4A*(T(i+1,4)-T(i,4)) - E4B*(T(i,4)-T(i,3)) - E4C*(T(i,4)-T(i,5));

    F(i,5) = -E5A*(T(i,5)-T(i-1,5)) + E5A*(T(i+1,5)-T(i,5)) - E5B*(T(i,5)-T(i,4)) + (-ep_5/(1-ep_5))*A56*(ep_5*sigma*(T(i,6)^4) +(1+ep_5)* ((ep_7*sigma*(T(i,6)^4)+(1+ep_7)*ep_5*sigma*(T(i,6)^4)) / (1-F22*(1-ep_7)*(1-ep_5)*F21)) -sigma*(T(i,5)^4));

    F(i,6) = -E7A*(T(i,6)-T(i-1,6)) + E7A*(T(i+1,6)-T(i,6)) + E7B*(T(i,6)-T(i,7)) + (-ep_7/(1-ep_7))*A67*(sigma*(T(i,6)^4) - (ep_7*sigma*(T(i,6)^4)+(1+ep_7)*ep_5*sigma*(T(i,5)^4)) / (1-F22*(1-ep_7)*(1-ep_5)*F21));

    F(i,7) = -E8A*(T(i,7)-T(i-1,7)) + E8A*(T(i+1,7)-T(i,7)) - E8B*(T(i,6)-T(i,7)) + E8C*(T_amb0-T(i,7)); 

    
       end
    end


    
end

Answer 1

[Shishir Reddy]

Hey Eoin

It looks like you're trying to model the temperature distribution in a superconducting cable using a finite volume method.

I see that you are accessing ‘T(i+1, : )’ when ‘i == N’, which would result in unexpected behaviour as index exceeds the array bounds. So, the boundary condition should be handled carefully at the last node (‘i == N’) as follows

elseif i == N
    % Boundary condition at the end of the cable
    F(i, 1) = mdot_int*cp_int*(T(i-1, 1) - T(i, 1)); 
    

Finally, ensure that initial conditions are set correctly and consistently for all nodes.

These changes should help address issues related to temperature jumps and ensure the stability of your numerical method.

I hope this helps.