plane wall with transient heat conduction and heat flux on one side that lasts for 0.2 seconds every 1 second, plot temperature distribution in matlab

To solve this problem, we need to use the finite difference method to discretize the governing equation and apply the temporal and spatial discretization techniques to calculate the temperature at each time step and position.

Assuming the wall is a homogeneous material and has a constant thermal conductivity, the governing equation for transient heat conduction can be written in the following form:

$$\rho C_p \frac{\partial T}{\partial t} = \frac{\partial}{\partial x}\left(k\frac{\partial T}{\partial x}\right) + q_x$$

where $\rho$ is the material density, $C_p$ is the specific heat capacity, $k$ is the thermal conductivity, $q_x$ is the heat flux on the wall surface, and $T$ is the temperature as a function of time and position.

We can discretize this equation using the central difference method to obtain:

$$\frac{T_{i}^{n+1}-T_i^{n}}{\Delta t}=\frac{k}{\rho C_p}\frac{T_{i+1}^n-2T_i^n+T_{i-1}^n}{\Delta x^2}+\frac{q_x}{\rho C_p}$$

where $n$ is the time index, $i$ is the spatial index, $\Delta t$ and $\Delta x$ are the time step and spatial step, respectively.

The heat flux $q_x$ is specified on one side of the wall and lasts for 0.2 seconds every 1 second. One way to implement this is to use a time-dependent boundary condition with a switch function:

$$q_x(t) = \begin{cases} q_0 & \text{if} , 0 \leq t \mod 1 < 0.2 \ 0 & \text{otherwise} \end{cases} $$

where $q_0$ is the constant heat flux value.

To solve this problem, we can follow these steps in Matlab:

1. Define the problem parameters: material properties, dimensions, initial and boundary conditions, time and spatial discretization.

main.m
% Material properties
rho = 1200; % density (kg/m3)
Cp = 1000; % specific heat capacity (J/kg-K)
k = 0.6; % thermal conductivity (W/m-K)

% Wall dimensions
L = 1; % length (m)
nx = 51; % number of spatial nodes
dx = L/(nx-1); % spatial step

% Initial temperature
T0 = 20*ones(nx,1); % uniform temperature (C)

% Boundary conditions
q0 = 5000; % heat flux (W/m2)
q = zeros(nx,1); % initialize boundary condition
tmax = 5; % maximum simulation time (s)
dt = 0.001; % time step

% Time-dependent boundary condition
q_switch = @(t) q0*(mod(t,1)<0.2);
548 chars
22 lines

2. Set up the coefficient matrix and the right-hand side vector for the finite difference equation.

main.m
% Coefficient matrix (tridiagonal)
A = zeros(nx);
A(1,1) = 1; % Dirichlet boundary condition (fixed temperature)
A(nx,nx) = 1; % Dirichlet boundary condition (fixed temperature)
for i = 2:nx-1
    A(i,i-1) = -k/dx^2;
    A(i,i) = rho*Cp/dt + 2*k/dx^2;
    A(i,i+1) = -k/dx^2;
end

% Right-hand side vector
b = zeros(nx,1);
323 chars
13 lines

3. Solve the equation iteratively using a time loop and the backslash operator.

main.m
% Time loop
for t = 0:dt:tmax
    % Update boundary condition
    q = q_switch(t);
    
    % Update right-hand side vector
    b(2:nx-1) = rho*Cp/dt*T0(2:nx-1) + q(2:nx-1);
    b(1) = T0(1); % Dirichlet boundary condition
    b(nx) = T0(nx); % Dirichlet boundary condition
    
    % Solve the system of equations
    T = A\b;
    
    % Update initial temperature for the next time step
    T0 = T;
end
405 chars
17 lines

4. Plot the temperature distribution at different times using the surf function.

main.m
% Spatial grid for plotting
x = linspace(0,L,nx);

% Time instants for plotting
tvals = [0 0.25 0.5 1 2 5];

% Temperature matrix for plotting
Tvals = zeros(nx,length(tvals));
for i = 1:length(tvals)
    [~,j] = min(abs(tvals(i)-linspace(0,tmax,length(T))));
    Tvals(:,i) = T(:,j);
end

% Plot temperature as a function of position and time
figure;
surf(x,tvals,Tvals');
xlabel('Position (m)');
ylabel('Time (s)');
zlabel('Temperature (C)');
444 chars
20 lines

This code will generate a 3D surface plot of the temperature distribution in the wall as a function of position and time. The plot will show how the temperature changes over time due to the periodic heat flux on one side of the wall.

Note that this code assumes a 1D problem with a uniform cross-section. To model more complex geometries, boundary conditions, or material properties, the code will need to be adapted accordingly.