The problem involves solving a 1-D heat transfer equation with a moving heat source using the finite difference method (FDM) with implicit formulation. The equation is:

$\frac{\partial (\rho c T)}{\partial t} + \frac{\partial (\rho c u T)}{\partial x} = \frac{\partial}{\partial x} \left(k \frac{\partial T}{\partial x}\right) + Q,$ where $Q(x, t) = 85e^{-2\frac{(x-ut)^2}{0.25^2}}$, and the boundary and initial conditions are given.

Step-by-step approach:

1. Discretize the domain:

   • Spatial domain: $x \in [0, 10]$.
   • Time domain: $t \in [0, 3]$.
   • Choose suitable grid spacing $\Delta x$ and $\Delta t$.

2. Apply finite difference approximations:

   • Use implicit scheme for time derivative.
   • Backward difference for the convective term $(\rho c u \frac{\partial T}{\partial x})$.
   • Central difference for spatial derivative of conduction term.

3. Formulate the system of equations:

   • Incorporate the boundary conditions:
     • $\frac{\partial T}{\partial x}(0, t) = 0$ translates to a Neumann boundary condition.
     • $T(10, t) = 20 + (e^t - 1)$.
   • Initial condition: $T(x, 0) = 20$.

4. Solve the linear system at each time step:

   • Use a matrix equation $A \mathbf{T}^{n+1} = \mathbf{b}$ where $A$ represents the coefficients from the discretization, $\mathbf{T}^{n+1}$ is the temperature at the next time step, and $\mathbf{b}$ includes known terms and boundary conditions.

Here is the MATLAB code for solving this problem:

matlabCopy code
% Parameters
c = 1.5e6; rho = 1e-6; k = 0.25; u = 0.8;
L = 10; t_max = 3;
Nx = 100; Nt = 300; % Number of grid points
dx = L / (Nx - 1); dt = t_max / Nt;
alpha = k / (rho * c);

% Discretized spatial and temporal grids
x = linspace(0, L, Nx);
t = linspace(0, t_max, Nt);

% Initialize temperature
T = 20 * ones(Nx, 1); % Initial condition
T_new = T;

% Implicit matrix coefficients
A = zeros(Nx, Nx);
b = zeros(Nx, 1);

% Fill the matrix A
for i = 2:Nx-1
    A(i, i-1) = -alpha / dx^2 - u / (2 * dx);
    A(i, i) = 1 / dt + 2 * alpha / dx^2;
    A(i, i+1) = -alpha / dx^2 + u / (2 * dx);
end
% Neumann boundary at x=0
A(1, 1) = 1; A(1, 2) = -1;
% Dirichlet boundary at x=L
A(end, end) = 1;

% Time stepping
for n = 1:Nt
    % Source term
    Q = 85 * exp(-2 * ((x - u * t(n)).^2) / 0.25^2);
    % Right-hand side
    b(2:end-1) = T(2:end-1) / dt + Q(2:end-1)';
    b(1) = 0; % Neumann BC at x=0
    b(end) = 20 + (exp(t(n)) - 1); % Dirichlet BC at x=L
    
    % Solve the system
    T_new = A \ b;
    
    % Update temperature
    T = T_new;
end

% Plot the results
[X, T_grid] = meshgrid(x, t);
surf(X, T_grid, T', 'EdgeColor', 'none');
xlabel('x'); ylabel('t'); zlabel('T(x,t)');
title('Temperature Distribution');

Notes:

1. Adjust $Nx$ and $Nt$ for desired spatial and temporal resolution.
2. Ensure stability and convergence by verifying $\Delta t$ and $\Delta x$ relationship.

Would you like a breakdown of any part of the code or further details on the methodology?

Questions to Expand:

1. How can we validate this numerical solution?
2. What are the key challenges in applying finite difference methods?
3. How do different boundary conditions impact the solution?
4. Could this problem be extended to 2-D or 3-D cases?
5. What are alternative numerical methods for solving such PDEs?

Tip:

Always check the stability criterion ($CFL$ condition) when using finite difference methods to ensure accuracy.