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

5 views (last 30 days)
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

Answers (1)

Shishir Reddy
Shishir Reddy on 13 Sep 2024
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.

Categories

Find more on Thermal Analysis in Help Center and File Exchange

Products


Release

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!