Order a Similar Paper Order a Different Paper
Write a function to solve the two-dimensional in Matlab,unsteady heat conduction equation with no internal heat generationon a square domain. The length of each side of the square is 1 m.The function will have the following 4 inputs: npts number of grid points in each coordinate direction nt number of time steps to take dt size of time step (s) alpha thermal diffusivity (m2/s) Use the following initial and boundary conditions: Initialize the body to T = 0 oC. At time t = 0,boundary conditions are imposed such that the temperatures on thesides of the square vary linearly. Here are the cornertemperatures: bottom left corner: T= 50 oC bottom right corner: T = 500oC top rightcorner: T = 350 oC top left corner: T = 150 oC Write the function definition statement. If you type
docfunction at the MATLAB command prompt it will give youdetails on the proper syntax. The syntax is general, showing how touse both inputs and outputs. We will have no outputs, only the fourinputs noted above. Make the name of your function the same as the name of yourfile. For example, if you call your function
heat,then your filename should be
heat.m. Create a variable name for each of the four corner boundaryconditions and the initial condition, and set their values. Noting that the temperature on each edge varies linearly,calculate the temperature at every grid point on the outerboundaries. You’ll need the corner temperatures to do thiscalculation. Create an array to contain the “old” temperature field at timelevel p, and call it told. A convenient way to do this is to usethe
ones function (look at doc ones), whichcreates an array of all ones of specified dimensions. An easy wayto set the field to the initial condition is to multiply this arrayby the initial temperature. If there are any calculated variables that you will use multipletimes, calculate them. For example, it’s handy to have a variablenpm = npts – 1, where npts is the input number of points in eachdirection. Create an array to contain the “current” solution at time levelp+1, and call it tnew. An easy way to do initialize it in MATLAB isto simply use tnew = told; Create a time loop going from 1 to nt, where nt is the inputnumber of time steps. Inside the time loop, create two other nested loops, one for the“m” indices (x direction) and one for the “n” indices (ydirection). Since we don’t have to re-calculate any values onboundaries, the loops will go from 2 to npm. For each pair of indices within the loops, calculate thesolution at time level p+1 (tnew) as a function of temperatures attime level p (told) using the finite-difference form of the 2-Dheat equation. We are going to do a simple animation of the unsteady results asthey evolve. Use the contourf function to create the evolvingcontour plot. (doc contourf) One of the arguments is the number ofcontours, and 20 seems to work well for this problem. We want tore-draw the contour plot once per time step, so contourf should becalled after the end of the nested spatial loops, but before theend of the time loop. It’s not necessary to use a plot handle as wedid in the example in class. Simply call contour each time throughthe loop. MATLAB tries to be efficient, so when it sees the contourffunction, it will buffer the results rather than plotting themright away. We don’t want it to buffer since we want to create areal-time animation. To prevent buffering, include a drawnowstatement on a separate line immediately after the call tocontourf. (This is similar to the example we did in class.) Now that the current time step is complete, we have to copy thenew solution to the old one so we can go on to the next iteration,told = tnew. Make this the last line before the end of the timeloop. Add the following to the end of the time loop. The first lineadds a “colorbar,” which is a scale showing the correspondencebetween temperatures and colors. The second line adds a label tothe colorbar. Finally, the third line gets rid of the plot axes,since we don’t need them here. hc = colorbar; hc.Label.String = ‘Temperature, deg. C’; axis off When you’re got your program ready to run, use the followinginput values: npts = 40 nt = 1000 dt = 0.1 alpha = 0.001 To run the program, simply type heat(40, 1000, 0.1, 0.001) . . .