ELEC0030 Numerical Methods 2019 Coursework Problem 2 Exercise 1 (10%) Consider a finite difference grid of regular spacing x∆ such that the boundary of the region of interest lies half way between grid lines. If c is a grid point next to the boundary, (as for example point 171 in the figure below), and a is one just outside of the boundary, derive an expression for the value of the function, u(x) at point a to use in the equation for the point c of the grid (for which a is a neighbour), when the condition that must be satisfied at the boundary is: ( )b b u h u u x ∞ ∂ = − − ∂ where h is a constant, b indicates the boundary and u∞ is a given value. Consider that in this case the function varies monotonically between point c and a (there is not a maximum or a minimum between them, that is, its derivative doesn’t change sign there). Exercise 2 (75%) Heat transfer in a solid body can be represented by the equation: 2u u t κ∂ = ∇ ∂ (1) with some appropriate boundary conditions where u is the temperature and κ is the thermal diffusivity. The figure shows the xy-cross section of a solid material body of a large dimension along the z axis with a regular finite differences grid superimposed. The material parameters and other conditions of the problem are considered as not varying in the z-direction so the problem can be approximated as a 2D problem. On the bottom surface in the figure, a temperature T0 is applied a t = 0 and is kept constant at all times after by a powerful heat source. The side and top boundaries shown with a thick black line are insulated and the insulation can be considered perfect. The surface portion marked by a dotted line is open to the air outside where there will be heat loss by convection1. The distant bulk average temperature in the air is Ta. Before the heat source is applied at the bottom, the temperature of the whole body is uniform and equal to Ta. Use finite differences and the mesh specified in the figure to find the temperature distribution as a function of time over the structure. 1 There will also be loss by radiation but we’ll leave that out in this problem. Finite differences grid ELEC0030 Numerical Methods 2019 Notes: The points on the bottom surface have not been numbered since their temperature values are known. A surface with a perfect insulation means that no heat flux can cross the surface. Heat flux density is ˆ uu n n ∂ −∇ ⋅ = − ∂ , where nˆ is the normal to the surface. Then, perfect insulation means: . . 0 i b u n ∂ = ∂ at all points on the insulated boundaries (i.b.). Convection is the heat transfer due to flow of a fluid. In this case, the air in contact with the hot surface. The magnitude of this flow will depend of the temperature difference between the surface and the bulk of the air, that is: [( ( , ) ]B B uk h u x y u n ∞ ∂ − = − ∂ (2) where h is the convective heat transfer coefficient, k is the thermal conductivity and u∞ is the average (bulk) air temperature at some distance from the contact surface. The differential equation (1) above can be normalised to have a non-dimensional function u(x,y,t) where the constant κ doesn’t appear explicitly ( 2u u t ∂ = ∇ ∂ ). Similarly, consider the normalised boundary condition (2) where h and k in are also one. Write a Matlab program that will do all the calculations and draw the plots. The frontal numbering scheme shown in the figure is convenient for the solution of the matrix problem. However, in a problem like this, it is probably more convenient to use row and column numbers (i and j) referring to the grid during the assembly of the matrices. Time stepping Although the simplest form to discretise the time derivative would be a forward difference, it will lead to instabilities. Consider instead using a Crank-Nicolson approach. From the discretised equation, you will notice that the matrices depend on the grid spacing and time step. If these are kept constant, the matrices will not change during the time stepping process and their assembly can be left out of the time stepping loop. However, short steps are needed at the start. If these are kept constant, it will mean a very large number of iterations and a long running time. Consider instead using a variable time stepping approach, gradually increasing the time step through the calculations. For example, the use of a grid spacing d = 0.05 and a starting time step of dt = 0.002 is suggested, with an increase of 15% per iteration: dt=dt*1.15. Stopping the iterations A tolerance of 810− on the absolute error is suggested to stop the iterations, when the error is defined as the norm of the difference between the solution vector u in two successive time steps: 1n nu u −− . Constructing the generic equation Note that the grid has 20 rows and 10 columns of nodes with unknown temperature values. The nodes are arranged in a regular rectangular pattern so their number (k) can be written in terms of their row and column numbers (i) and (j). (Warning: Do not confuse these grid row and column numbers with those of the resultant matrices! – the order the matrices will be 200: 200 rows and 200 columns). Find an expression to convert between the grid row-column notation to the global numbering easily. Assign the row numbers to the grid of numbered points from the bottom (row 1: points 1 – 10, and so on and for the columns, from left to right (col. 1 contain the points 1, 11, 21, …). ELEC0030 Numerical Methods 2019 Assembly of the matrices: Write down in you report the general form of the algebraic equation that results from applying the differential equation (1) to a generic node in the structure, for example the one corresponding to node 146, that will generate the values for row 146 (of length 200) of the matrices A and B. Identify all nodes for which this equation applies (internal nodes) by their row and column numbers in the grid and implement that in your program to fill the matrices. Consider then all special cases one by one: Nodes next to the Dirichlet boundary (2 – 9), those next to the side boundaries (11, 21, … 121, 20, 30, 40, …190), those next to the top boundary (192 – 199). For convenience, treat the corner nodes (1, 10, 191, 200 – where two boundary conditions apply), separately. List these different cases in your report giving the corresponding modified equation. You will notice that when deriving the discretised form of the equation for all these cases a matrix equation can be written of the form: 1n n+ = +Au Bu V Write a piece of code to fill the corresponding values in each case, separating each segment by comment lines indicating what is intended. Notice that the matrices are sparse and have the same sparsity pattern. Sparse storage is then strongly recommended (top marks will only be available if this route is followed). In this case, you will have to fill arrays for the row and column indices and for the corresponding matrix values. (Look for the Matlab documentation for the command “sparse”). The vector V will still be defined as full and needs to be initialised with: V=zeros(200,1). If you decide to assemble full matrices instead (but you won’t get full marks), start by pre-allocating space for the matrices A and B as: A=zeros(200,200) to overwrite later with the values that are not zero. Once the matrices and the vector are assembled, use the command spy(A) to examine the sparsity pattern of the matrices. It should be the same for both matrices, verify this and copy this figure in your report. Does this pattern agree with the expected shape, given the expressions used to calculate the temperature? Comment on this, referring also to the total number of places in the matrix and the number of nonzero values. Run the program with the values: Time step: ∆t = 0.002 (initial value) Grid spacing (equal for x and y) d: =
0.05 Initial temperature and bulk temperature of the air: Ta = 10° Heat source temperature: T0 = 80°. Output of results and plots Apart from the sparsity pattern, plot the variation of temperature versus time as the iterations progress. For this, select a sample of points: For example: 5, 15, 25, … That is, those in column 5, augmented with the boundary values. That will give a sequence of line plots showing the temperature profile along that line during the iterations. Label the axes and put a title. At the end of the iterations, draw a surface plot of the final temperature distribution. For this, you will have to create a 2D array from the results vector u using the Matlab command “reshape” and augment that array with extra rows and columns to include the values at the boundaries. Produce also a contour plot of the final temperature distribution over the cross-section of the structure. Exercise 3 (15%) Consider the same situation described in Exercise 2, but now, we are only interested in a steady state solution. How would eqn. 1 change in this situation? (Hint: The clue is in the name: ELEC0030 Numerical Methods 2019 steady state). Write down the discretised equation for a generic node in this case and explain how all the other equations would change. Modify a copy of the program written for Exercise 2 to solve this case. Compare the solution you get now with the plots of the results found at the end of the time stepping in Exercise 2. Submit your report through the Moodle page as a pdf file. Write your name in the front page. The report should be a detailed account of the procedure chosen for these calculations, the results obtained and the comments or discussion requested. It should include the plots asked and the programs used for the calculations. You may refer to and cite parts of the code in the text if that is necessary to explain something but submit the complete program for each exercise as an appendix to the report. Do not include Introduction, summary, conclusions or theory sections; just report what you have done including the equations and necessary derivations and include as appendices the programs used for the three exercises. Marks will be allocated to the correctness and completeness of the answers, the presentation of your report and the quality and clarity of the programs written.
