create a temperature distribution for a plane wall under transient heat conduction in matlab

main.m
% Define the geometry and initial conditions of the wall
L = 0.1;     % length of wall in meters
W = 0.1;     % width of wall in meters
D = 0.01;    % thickness of wall in meters
k = 10;      % thermal conductivity of wall in W/(m*K)
c = 1000;    % specific heat of wall in J/(kg*K)
rho = 1000;  % density of wall in kg/m^3
T0 = 20;     % initial temperature of wall in degC
T1 = 100;    % temperature on left boundary in degC
T2 = 20;     % temperature on right boundary in degC

% Discretize the wall into small elements
nx = 50;  % number of elements along x direction
ny = 50;  % number of elements along y direction
dx = L/nx;
dy = W/ny;
[x,y] = meshgrid(dx/2:dx:L-dx/2, dy/2:dy:W-dy/2);

% Define time parameters
tmax = 1000;       % maximum simulation time in seconds
dt = 0.1;          % time step size in seconds
nt = round(tmax/dt); % number of time steps

% Initialize temperature matrix
T = T0*ones(nx,ny);

% Apply finite difference method to solve heat conduction equation
for n = 1:nt
    % Compute spatial derivatives using central difference scheme
    Tx = (T(2:end,:)-T(1:end-1,:))/dx;  % dT/dx
    Ty = (T(:,2:end)-T(:,1:end-1))/dy;  % dT/dy
    Txx = (T(3:end,:)-2*T(2:end-1,:)+T(1:end-2,:))/dx^2;  % d^2T/dx^2
    Tyy = (T(:,3:end)-2*T(:,2:end-1)+T(:,1:end-2))/dy^2;  % d^2T/dy^2
    
    % Compute temperature at next time step using backward Euler method
    Tnew = T + dt*k/(rho*c)*(Txx + Tyy);
    Tnew(:,1) = T1;   % apply left boundary condition
    Tnew(:,end) = T2; % apply right boundary condition
    T = Tnew;
end

% Plot the temperature distribution at final time step
surf(x,y,T')
xlabel('x (m)')
ylabel('y (m)')
zlabel('Temperature (degC)')
1677 chars
47 lines