I'm terrible with Matlab and am trying to find a nodal temperature distribution of a shape and cannot get my code even close to doing what I need

Question

Lee Brunton 2023-5-9

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1960184-i-m-terrible-with-matlab-and-am-trying-to-find-a-nodal-temperature-distribution-of-a-shape-and-canno

评论： Walter Roberson 2023-5-9

clear;
tic
k = 5; %Thermal conductivity of concrete - W/mK
T_w = 5; %Temperature of water - DegC
%Temperature surface 
h = 100; %Convective heat transfer of water - W/m^2K
q_gen = 500; %Heat generation from top surface W/m^2
%T_i = 40; %Intial guess of concrete temperature - DegC
R = 5; %Grid resolution - mm
dx = 0.001*R; %Distance between nodes in x direction - m
dy = 0.001*R; %Distance between nodes in y direction - m
A = (100/R)+1;
B = (60/R)+1;
C = (30/R)+1;
D = (20/R)+1;
E = 5/R+1;
F = 35/R+1;
G = 50/R+1;
H = 25/R+1;
I = 75/R+1;
AA = 80/R +1;
T = ones(B,A)*T_w;
Tp = T;
%x = 1;
i = 0;
error = 1;
max_error = 1e-6;
max_iter = 10000;
q = zeros(1,A);
Ts = zeros(1,A);
L = (0:R:100);
N = 1;
m = 1;
n = 1;
Temp = [];
while error>max_error && i<max_iter
    i = i + 1;
    Tp = T;
    for m = A:1
        for n = 1:B
            
            %Water temperature ("Node 13")
            if ((m==H)&&(m==I))&&((n==F)&&(n<G))
                if i==1
                    Tp(m,n) = T_w;
                end
                %Top left corner Node ("Node 1")
            elseif (m==1)&&(n==1)
                T(m,n) = (Tp(m,n+1)+Tp(m,n-1)+Tp(m+1,n))/3 + (2*q_gen*dx*dx)/(3*k);
                
                
                %Top horizontal edge of concrete ("Node 2")
            elseif ((1==m)&&(m<B))&&(n==1)
                
                
                T(m,n) = (Tp(m-1,n))/6 + (Tp(m,n+1)+Tp(m,n-1)+Tp(m+1,n))/3 + (q_gen*dx^2)/(3*k);
                
                
                %Fully Enclosed Concrete Nodes ("Node ")
            elseif (1<m)&&(m<B)&&(1<n)&&(n<D)|| (1<m)&&(m<B)&&(1<n)&&(n<B)
                T(m,n) = (Tp(m-1,n)+ Tp(m+1,n)+Tp(m,n+1)+Tp(m,n-1))/4;
                
                %Left vertical edge  %("Node ")
            elseif (m==1)&&(1<n)&&(n<A)
                T(m,n) = (Tp(m,n+1)+Tp(m,n-1))/4 + (Tp(m+1,n))/2;
                
                %Top left Corner of Water Node 59
            elseif (m==D)&&(n==C)
                
                T(m,n)= (Tp(m-1,n)+Tp(m,n+1))/(3+(h*dx)/k)+ (Tp(m+1,n)+Tp(m,n-1))/(6+(2*h*dx)/k) + (h*dx*T_w)/(3*k+h*dx);
                
                %Top Right Corner of Water
                %
            elseif (m==D)&&(n==G)
                
                T(m,n) = (Tp(m+1,n)+Tp(m,n+1))/(3+(h*dx)/k)+ (Tp(m-1,n)+Tp(m,n-1))/(6+(2*h*dx)/k) + (h*dx*T_w)/(3*k+h*dx);
                
                
                %Bottom Left Corner of Water
            elseif (m==AA)&&(n==C)
                
                T(m,n) = (Tp(m-1,n)+Tp(m,n-1))/(3+(h*dx)/k)+ (Tp(m+1,n)+Tp(m,n+1))/(6+(2*h*dx)/k) + (h*dx*T_w)/(3*k+h*dx);
                
                %Bottom Right Corner of Water Node 219
                
            elseif (m==AA)&&(n==G)
                T(m,n) = (Tp(m+1,n)+Tp(m,n-1))/(3+(h*dx)/k)+ (Tp(m-1,n)+Tp(m,n+1))/(6+(2*h*dx)/k) + (h*dx*T_w)/(3*k+h*dx);
                %Top of Water
            elseif (m==D)&&(C<n)&&(n<G)
                T(m,n) = (Tp(m,n+1)*k + (Tp(m-1,n)+Tp(m+1,n))*(k/2-(h*dx)/2)+h*dx*(Tp(m,n-1)+T_w))/(2*k+h*dx);
                %Bottom of Water
                
            elseif (m==AA)&&(C<n)&&(n<G)
                
                T(m,n) = (Tp(m,n-1)*k + (Tp(m-1,n)+Tp(m+1,n))*(k/2-(h*dx)/2)+h*dx*(Tp(m,n+1)+T_w))/(2*k+h*dx);
                %Left side of Water
            elseif (D<m)&&(m<AA)&&(n==C)
                T(m,n) = (Tp(m-1,n)*k + (Tp(m,n+1)+Tp(m,n-1))*(k/2-(h*dx)/2)+h*dx*(Tp(m+1,n)+T_w))/(2*k+h*dx);
                %Right Side of Water
                
            elseif (D<m)&&(m<AA)&&(n==G)
                T(m,n) = (Tp(m+1,n)*k + (Tp(m,n+1)+Tp(m,n-1))*(k/2-(h*dx)/2)+h*dx*(Tp(m-1,n)+T_w))/(2*k+h*dx);
                %Bottom Nodes Insulated
            elseif (m==A)&&(1<n)&&(n==B)
                T(m,n) = (Tp(m-1,n)+Tp(m+1,n))/4 + (Tp(m,n+1))/2;
                %Bottom Left Corner Node
            elseif (m==A)&&(n==1)
                T(m,n) = (Tp(m+1,n)+Tp(m,n+1))/2;
                
                
            end
        end
    end
    
    %Error Calculation
    error = max(abs(T-Tp),[],'all');
    it_error(i) = error;
    Tp = T;
end
T = T';
for N = 1:B
    q(N) = h*dx*(T(1,N)-T_w);
    
    Ts(N) = T(1,N);
end
q_total = sum(q);
Q = q_total/0.1;
%Plotting the Temperature Distributions
% figure(1)
% imagesc(T);
% colorbar
% title('Temperature Distribution of Section')
% figure(2)
% contourf(flipud(T), 10)
% colorbar
% title('Temperature Distribution of Section in 10 levels')
% figure(3)
% surfc(T')
% colorbar
% title('Topological Graph of Section')
% figure(4)
% plot(L,Ts)
% xlabel('Distance from top left corner (mm)')
% ylabel('Temperature (Degrees Celcius)')
% title('Temperature Profile Over Surface')
% figure(6)
% loglog(1:i, it_error)
% xlabel('Iteration Number')
% ylabel('Residual Error')
% title('5mm Mesh Error Convergence Plot')
% toc