solve liquid phase mass balance and energy on a spatially 2d domain in cylindrical coordinates in matlab

main.m
% Parameters
r_min = 0;
r_max = 1;
z_min = 0;
z_max = 2;
nr = 20;
nz = 40;
rho = 1000; % Density
c_p = 4182; % Heat capacity
k = 0.6; % Thermal conductivity
h = 10000; % Heat transfer coefficient
T_w = 373; % Wall temperature
T_i = 293; % Initial temperature

% Grid
dr = (r_max-r_min)/(nr-1);
dz = (z_max-z_min)/(nz-1);
r = linspace(r_min,r_max,nr);
z = linspace(z_min,z_max,nz);
[RR,ZZ] = meshgrid(r,z);
r_c = r(1:end-1)+dr/2;
z_c = z(1:end-1)+dz/2;
[RR_c,ZZ_c] = meshgrid(r_c,z_c);

% Time steps
t_end = 1000;
dt = 1;
t = 0:dt:t_end;
nt = length(t);

% Initial conditions
T = T_i*ones(nr,nz);

% Mass balance
for n = 1:nt-1
    % Mass fluxes
    J_r = -k*rho.*diff(T,1,1)./dr;
    J_z = -k*rho.*diff(T,1,2)./dz;
    
    % Mass balance
    dM_r = 2*pi*r_c(1)*dz*rho.*J_r(:,1);
    dM_z = 2*pi*r_c'*dr*rho.*J_z(1,:);
    dM_r(2:end) = dM_r(2:end) - dM_r(1:end-1);
    dM_z(:,2:end) = dM_z(:,2:end) - dM_z(:,1:end-1);
    M = cumsum(dM_r,1)+cumsum(dM_z,2);
    
    % Update temperature
    dT_r = diff(J_r,1,1)./dr;
    dT_z = diff(J_z,1,2)./dz;
    T_c = T(1:end-1,1:end-1)+.25*(dT_r(:,1:end-1)+dT_r(:,2:end)...
        +dT_z(1:end-1,:)+dT_z(2:end,:))*dt;
    T = T*0;
    T(2:end-1,2:end-1) = T_c;
    
    % Boundary conditions
    T(:,1) = T_w;
    T(1,:) = T(2,:);
    T(end,:) = T(end-1,:);
end

% Visualize 
figure
contourf(RR,ZZ,T)
colorbar
xlabel('r [m]')
ylabel('z [m]')
title('Temperature [K]')
1408 chars
68 lines