Lab 1 -- Solving a heat equation in Matlab
- To introduce students to programming, and Matlab programming in particular.
- To learn how pseudocode is used as an intermediate step in converting an algorithm stated in English and Math into a computer program.
- To operationalize calculus concepts covered in lecture.
- To introduce finite difference methods for the numerical solution of differential equations.
Consider the problem of determining the subsurface temperature fluctuations of rock as the result of daily or seasonal temperature variations. As the surface is heated or cooled, the heat will diffuse through the soil and rock. Mathematically, diffusion may be represented by the equation
where T(t,z) is the temperature at time t and depth z, and K is a constant which measures the "diffusivity" of the rock, i.e. how quickly heat moves through the rock. A typical value for K is 10^-6 m^2/s.
- Question 1: Interpret the heat equation in words. What is T(t,z)? dT/dt? dT/dz? d^2T/dz^2? Why should the left hand side and right side equate?
The above equation may be solved "analytically", i.e. T(t,z) may be written exactly as an equation. Matlab is more appropriate, however, for solving problems numerically, i.e. approximately as (a matrix or a map) of numbers.
In order to solve the system numerically, we must first choose an appropriate "space" to solve the problem in. We will have two axes, the vertical axes representing depth, from 0m to 100 m below the surface, and the horizontal axis will represent one year's worth of time.
The heat equation is an example of what is known as a "partial differential equation." A differential equation is any equation in which a function (Temperature in time and space in this instance) is not represented directly, but via it's derivative. Partial indicates that there are at least two variables (time and space) in the derivatives. An ordinary DE only has one derivative.
In order to solve a PDE numerically, we need to specify boundary conditions. The four boundaries of our "space" are the surface (0m), the bottom depth (100m), an arbitrary beginnning time, and one year later.
- Top Boundary:
The top boundary is the simplest, we will simply specify a temperature that varies seasonally, i.e.
This will create an average annual temperature of 15, an annual low of 5 and an annual high of 25 (time is measured in months, and we are not starting on Jan 1).
- Bottom Boundary:
The bottom boundary condition is easy as well... we will assume that no heat reaches the bottom from the top, i.e.
This specifies that there is 0 heat flow at the bottom. We do not yet know what this bottom temperature will be.
- Left and right boundaries
The left and right boundary conditions are a bit more problematic, as they require that we already know the temperature at every depth, which is exactly what we are trying to find out! But we do know that the temperature now and exactly one year from now should be the same.
The final step is to determine how to use finite differences to evaluate the derivatives. For the first derivative of temperature with respect to time we will use the "forward derivative"
For the first derivative of temperature with respect to depth, we will use the "central derivative"
% For the second derivative of temperature with respect to depth, we will use
See your class notes for the logic use here.
The following little program implements this procedure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Script to solve the heat equation in soil with %%% seasonal temperature forcing %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%% Initialize variables dz = .25; %each depth step is 1/4 meter Nz = 400; % Choose the number of depth steps (should go to at least 100 m) Nt = 5000; % Choose the number of time steps dt = (365*24*60*60)/Nt; %Length of each time step in seconds (~ 6.3*10^3 seconds, or ~105 minutes) K = 2*10^-6; % "canonical" K is 10e-6 m^2/s T = 15*ones(Nz+1,Nt+1); %Create temperature matrix with Nz+1 rows, and Nt+1 columns % Initial guess is that T is 15 everywhere. time = [0:12/Nt:12]; T(1,:) = 15-10*sin(2*pi*time/12); %Set surface temperature maxiter = 500 for iter = 1:maxiter Tlast = T; %Save the last guess T(:,1) = Tlast(:,end); %Initialize the temp at t=0 to the last temp for i=2:Nt+1, depth_2D = (T(1:end-2,i-1)-2*T(2:end-1,i-1)+T(3:end,i-1))/dz^2; time_1D = K*depth_2D; T(2:end-1,i) = time_1D*dt + T(2:end-1,i-1); T(end,i) = T(end-1,i); % Enforce bottom BC end err(iter) = max(abs(T(:)-Tlast(:))); %Find difference between last two solutions if err(iter)<1E-4 break; % Stop if solutions very similar, we have convergence end end if iter==maxiter; warning('Convergence not reached') end
maxiter = 500
- Question2: Write pseudocode for the above algorithm and turn it in.
figure(1) plot(log(err)), title('Convergence plot') figure(2) imagesc([0 12],[0 100],T); title('Temperature plot (imagesc)') colorbar figure(3) depth = [0:dz:Nz*dz]; contourf(time,-depth,T); title('Temperature plot (contourf)') colorbar
figure(4) plot(time,T(1,:),time,T(21,:),time,T(41,:),time,T(61,:),time,T(81,:)) xlabel('Time (months)'); ylabel('Temperature (C)'); legend('0m','5m','10m','15m', '20m')
- Question 3: Modify the above code to evaluate heat diffusion for daily forcing. Don't forget to change the time and depth steps to appropriate values. (You probably want to keep 5000 time steps and 400 depth steps, however. You may also want to change the amplitude of the temperature.
Plot the solution using contourf. Adjust your plot so that only the "interesting" depths are shown. Plot the solution at 0, 1/2, 1 and 1.5 meters below surface. About how many hours after the surface reaches its hottest, does it reach its hottest at 1/2 meters below the surface?