I'm trying to solve a heat diffusion problem (with finite difference approximation) with conditions stated like this:
A 10cm thin plate is initially at $120^{\circ}C$ then suddenly the right side is dropped to $0^{\circ}C$ while the left surface is kept insulated. Given $\alpha$, find the temperature distribution at $10$ and $20$ seconds.
I'm supposed to use the Crank-Nicolson approximation, but I'm not sure how to construct matrices $A$ and $b$ to solve $Ax=b$.
I'm pretty sure the "left surface insulated" means (in MATLAB notation) :
A(1,:) = [1, -1, zeros(1,Nx-2)];
b(1) = 0;
So for each time step, $T_1 - T_2 = 0$, i.e. there is no heat flow.
But I'm not sure how to interpret "right side suddenly dropped to 0". I did it like this:
Solution = zeros(Nt, Nx); % Time steps x Spatial nodes
Solution(1,:) = [ones(1,Nx-1)*120, 0]; % so the initial distribution is [120, 120, ... , 120, 0]
Is this correct? And does the question imply right surface is kept at $0^{\circ}C$? So I keep solving the indermediate steps as :
A(Nx,:) = [zeros(1,Nx-1), 1];
b(Nx) = 0;
i.e. at each time step $T_{(N_x)} = 0$.
Does this look correct? The problem is, I get an unstable solution, whereas CN method is said to be "unconditionally stable". The solution oscillates, and also nodes very close to right hand side become $<0$ so there is a sharp jump at the right hand side. See:
What is surprising is, as I increase Nx, therefore impose a finer grid, the solution gets worse, opposite of what I expected. This screenshot is from my trial with 100 points, and T_right = 10 instead of 0.
Thanks for any help!
As long as I understand you should not have a jump. What is your grid step? In general the finer the grid, the better the result, but numerically you cannot refine it forever.
Assume the problem on an interval $[a,b]$.
The initial condition is that all plate has $120^\circ$ at $t_0$.
Insulated boundary condition for heat diffusion problem means no flow heat flow, i.e. $u_x(a,t)=0$.
The other side is a different story, you have to use the Newton's Law of Cooling (see also), which (for the second end) says
$$-K_0 u_x(b,t)=H (u(b,t) - u(b,t_0))$$ where $K_0$ is thermal conductivity (in $u_t = k u_{xx}$ you have thermal diffusion $k=\frac{K_0}{c\rho}$, $c$- specific heat, $\rho$- mass density) and $H>0$ is some coefficient.
In your case it looks like, you should use the following BC
$$u_x(b,t)=\beta (u(b,t) - 120^\circ)$$ for some $\beta$ (not sure if the $\alpha$ you talking about is $k$ or this $\beta$).
For implementation example look
http://web.cecs.pdx.edu/~gerry/class/ME448/notes/pdf/Alt_BC_slides.pdf