MEC3456 Assignment Page 1 of 9 MAE3456/MEC3456 LAB 06 Due: 11:55PM, Wednesday 26th May 2021 (Middle of Week 12) This lab should be completed INDIVIDUALLY. Plagiarism will result in a mark of zero. Plagiarism includes letting others copy your work and using code without citing the source. Collaborating with others to discuss algorithms and details of MATLAB syntax and structures is acceptable (indeed encouraged), however you MUST write your own MATLAB code. All assignments will be checked using plagiarism-detecting software and similarities in submitted code will result in a human making a decision on whether the similarity constitutes plagiarism. INSTRUCTIONS Download template.zip from Moodle and update the M-Files named Lab06_Q1.m, etc… with your Lab code. DO NOT rename the M-Files in the template or modify run_all.m Check your solutions to the questions by running run_all.m and ensuring all questions are answered as required. SUBMITTING YOUR ASSIGNMENT Submit your assignment online using Moodle. You must include the following attachments: A ZIP file (NOT .rar or any other format) named in the following way: Surname_StudentID_Lab_6.zip (e.g. Rudman_23456789_Lab_6.zip) The zip file should contain one folder containing everything needed for Lab06, i.e. the following files: a. All MATLAB m-files for lab tasks: run_all.m, LabNN_ Q1.m, etc… b. Any additional MATLAB function files required by your code c. All data files needed to run the code, including any input data provided to you d. Any hand calculations or written responses asked for - scanned in as a SINGLE PDF file YOUR ZIP FILE WILL BE DOWNLOADED FROM MOODLE AND ONLY THOSE FILES INCLUDED IN YOUR SUBMISSION WILL BE MARKED – YOU WILL GET ZERO MARKS FOR ANY MISSING FILES – DON’T ASK ME TO ACCEPT THEM LATE AS I WILL SAY NO. We will extract (unzip) your ZIP file and mark the lab based on the output of run_all.m and any hand calculations or written responses. It is your responsibility to ensure that everything needed to run your solution is included in your ZIP file. STRONG RECOMMENDATION: After uploading to Moodle, download your Lab to a NEW folder, extract it and open a NEW instance of MATLAB. Ensure your code runs as expected. CHECK YOUR PDF IS INCLUDED. You can upload your submission as many times as you like BEFORE the submission deadline. After this time, it is not possible to resubmit your lab without incurring a late penalty. If you forget to include some of your code or your PDF you cannot include it later without penalty. MEC3456 Assignment Page 2 of 9 MARKING SCHEME This lab is marked out of 75 and full marks is worth 8% of your total unit grade for the semester. This lab will be graded using the following criteria: 1) run_all.m produces results automatically (no additional user interaction needed except where asked explicitly – NOTE, I have included pause commands in run_all so that intermediate answers can easily be viewed by the demonstrators – please don’t remove them) 2) Your code produces correct results (printed values, plots, etc…) and is well written. 3) Your answers to the hand-written questions are correct 4) Programming style, efficiency of algorithm and quality of output (figures, tables, written text ...) will be assessed in this lab and marks will be lost for poor style. ASSIGNMENT HELP 1) You can ask questions in the Discussion Forum on Moodle 2) Hints and additional instructions are provided as comments in the assignment template M-Files 3) Hints may also be provided during lectures 4) The questions have been split into sub-questions. It is important to understand how each sub- question contributes to the whole, but each sub-question is effectively a stand-alone task that does part of the problem. Each can be tackled individually. 5) I recommend you break down each sub-question into smaller parts too, and figure out what needs to be done step-by-step. Then you can begin to put things together again to complete the whole. 6) To make it clear what must be provided as part of the solution, I have used bold italics and a statement that (usually) starts with a verb (e.g. Write a function ..., Print the value..., etc.) SUBMISSION CHECKLIST Check off these boxes AFTER you have submitted your Lab 1) I have downloaded my Lab from Moodle and unzipped it into a NEW FOLDER. o 2) I have checked the PDF of hand calculations is included, and if missing I have added it. o 3) I have checked my code runs as required and fixed any errors such as missing files, etc. o 4) If any changes to my submission were required from steps 2 or 3 above, I have uploaded my submission AGAIN and undertaken steps 1-3 AGAIN. o MEC3456 Assignment Page 3 of 9 QUESTION 1 [25 MARKS TOTAL] Background In this question we consider the 1-D advection-diffusion equation: Eqn (1) In this equation, c is a “velocity” (assume that c > 0) and k is a diffusion coefficient (k > 0). You will solve this equation on a domain 0 ≤ x ≤ L with boundary conditions at x = 0 given by q(x=0,t) = 0 and ∂q/∂x(x=L,t) = 0. The initial condition is a square wave given by Q1a By hand, write the forward in time, centred in space (FTCS) difference equation that approximates equation (3). What is the temporal and spatial order of accuracy of this discretization? Q1b Undertake a von Nuemann stability analysis for this difference equation (you can assume that the error satisfies the same difference equation as the solution for q, because it does). Consider just ONE Fourier mode of the representation of the error, and write it in the simplified form Write an expression for the amplification factor (i.e. the ratio of error at time n+1 to that at time n) and replace any complex exponential terms by expressions involving sin and cos - it will be a complex number. Simplify as far as you can, but you do NOT need to derive the stability criterion. Q1c Undertake a Taylor series stability analysis - i.e. expand each term in the finite difference equation about time n, location j, as done in workshops and replace the time derivative operator with a spatial derivative operator that looks like . Then examine what the coefficient of the is in the resulting PDE and determine the timestep restriction for stability. ∂q ∂t + c ∂q ∂x = k ∂2q ∂x2 q(x, 0) = 1 if 15< x < 200 otherwise ! " # xiKjnn j eA D=e ∂ ∂t = −c ∂ ∂x + k ∂ 2 ∂x2 ∂2q ∂x2 MEC3456 Assignment Page 4 of 9 Q1d Write a MATLAB function that updates the values of q by ONE TIMESTEP using the FTCS scheme (it will be a short function). Ensure that your function has the function header function [q_new] = ftcs_step(q,dt,dx,c,k) The boundary condition at x=0 is easily implemented using q(x=0)=0, but you will need to set the boundary condition at x=L using a finite difference approximation to the derivative. Thus, when q(L,t) is needed in an expression, you will need to use the finite difference expression of the derivative to set the BC. (In this lab, centered or one-sided differences are OK to use for the BC - whatever you like). Q1e Modify the template in Lab06_Q1.m using the following values • Set L = 100 • Discretise your domain with 250 intervals. • Choose c = 1.0 • Set tmax = 60 Choose 3 values of k [0.01, 0.1, 1.0] and integrate forward in time until t = tmax with the following CFL1 numbers: [0.01, 0.1, 0.25] and (i.e. you will have 9 solutions). Create a single plot with 9 subplots. Each subplot will have the solution for one CFL number and one value of k. In each subplot, plot the solution at approximately every 20 time units (i.e. at t = 0, 20, 40, 60) on the same subplot (i.e. you should have 9 subplots, each with 4 curves). The plot times can be “approximate” (e.g. 19.95 is close enough to 20). ALSO include a plot of the analytic solution of the advection-only problem at each time in a different colour in each subplot for comparison (giving 8 curves per subplot). Describe the key features of what you are seeing as you change k and the CFL number and suggest WHY you observe what you observe (the following hints might help you). HINTS: 1) The FTCS scheme for equation 3 when k=0 (i.e. pure advection) is unconditionally unstable. 2) The FTCS scheme for equation 3 when c = 0 (i.e. pure diffusion) is conditionally stable. The condition for stability is discussed in the lecture notes. 3) The answer to Q1c could provide some insight. 1 The CFL (or Courant) number is defined as (c ∆t/∆x) and is related to what fraction of a grid cell “information” is advected in one timestep. MEC3456 Assignment Page 5 of 9 QUESTION 2 [15 MARKS TOTAL] Background The transient heat equation in cylindrical coordinates (, , ) may be described using: Eqn (2) where is temperature and > 0 is the thermal diffusivity. This is generally a 3D problem. However, for the case when the domain and boundary conditions are axisymmetric, Eqn (2) reduces to the 2D cylindrical coordinate representation given in Eqn 3: Eqn (3) Q2a Write the BTCS (implicit) scheme for the 2D axisymmetric heat equation in Eqn 3 at node (j,k) for time “n+1”. The index j represents location in the r-direction, while the index k represents location in the z-direction. You must move any unknowns (at time “n+1”) to the LHS and all known values (at time “n”) to the RHS. Q2b Undertake a von Neumann stability analysis of the difference equation for the BTCS scheme in Q2a and derive an expression for the magnitude of the ratio of errors at time n+1 to those at time n. You should write the error as where R and S are in the r- and z-directions, respectively. You can treat R and S in the same way we treat K (=2pk/L) in the 1D heat equation. Is the BTCS scheme is the scheme unstable, conditionally stable or unconditionally stable? Justify your answer. If it is conditionally stable, derive the timestep restriction for stability. Q2c We can write the BTCS scheme in Q2a as a matrix problem of the form where Q is the temperature vector, [A], [B] are matrices that depend on the identity matrix [I], the Laplacian matrix [L], the timestep Dt and the diffusion coefficient a. BC is a vector that stores values associated with the boundary conditions. Specify exactly what the entries of [L] are equal to in terms of r, dr, dz and what [A], [B] are equal to in terms of [I], [L], Dt, and a (you can ignore boundary conditions in your answer). ∂Q ∂t =α 1 r ∂ ∂r r ∂Q ∂r ⎛ ⎝⎜ ⎞ ⎠⎟ + 1 r 2 ∂2Q ∂θ 2 + ∂ 2Q ∂z2 ⎧ ⎨ ⎩ ⎫ ⎬ ⎭ ∂Q ∂t =α 1 r ∂ ∂r r ∂Q ∂r ⎛ ⎝⎜ ⎞ ⎠⎟ + ∂ 2Q ∂z2 ⎧ ⎨ ⎩ ⎫ ⎬ ⎭ ε jk = A neiRjδ reiSkδ z Qn+1 = [A]−1 [B]Qn +BC( ) MEC3456 Assignment Page 6 of 9 QUESTION 3 [25 MARKS TOTAL] Consider a flat hollow cylinder of inner radius cm and outer radius 2 cm. We will consider a case which is axisymmetric (i.e. does not depend on q) and is shown in Figure 1. Figure 1: A hollow cylinder (left) modelled as a 2D axisymmetric geometry (right) The temperature distribution (, , ) across the cylinder may be described using the axisymmetric heat equation given in Eqn (3). The cylinder has the following boundary condition: (, , ) = sin(2) Eqn (4) and initial condition: (, , 0) = sin(2) Eqn (5) Q3a In your PDF, outline the algorithm needed to solve Eqn 3 using the BTCS scheme from time = 0 to time = tn, with a constant timestep, Dt (assume that tn = NDt). Include the major steps required to set the problem up (don’t include minor steps like setting grid spacing, etc.) This algorithm must be in the form of pseudo code or a flow chart. An example of pseudo code for solving an ODE using the Euler method is shown in the snippet below. NOTE: All BCs in this problem are Dirichlet conditions (i.e. QBC is known in advance). THUS, it is possible to remove the Boundary points from the set of unknowns, which is what is recommended below. If you want to keep them as unknowns, then the size of the grid will be Nr*Nz instead of (Nr-2)*(Nz-2) as specified in Q3b. Set up parameters of the problem Determine stepsize and number of steps for each time step Solve the Euler method Store calculated values in vector end MEC3456 Assignment Page 7 of 9 Q3b Write a Matlab function that sets up the Laplacian matrix L and the boundary condition vector BC. Your function header must be written as function [L,BC] = MatSet(Nr,Nz,rg,zg,alpha,BCf) where • Nr and Nz are the number of points in the r- and z-directions, respectively • rg and zg are column vectors containing the r- and z-values of the grid points respectively. • a is the thermal diffusion coefficient • BCf is a function describing the boundary condition as a function of r,z and t (see Eqn 4) • L is the ((Nr-2)*(Nz-2) x (Nr-2)*(Nz-2)) Laplacian matrix • BC is a ((Nr-2)*(Nz-2) x 1) column vector on the RHS that accounts for the boundary conditions NOTE: You can define an anonymous function that converts a 2D node number (j,k) = (jdr, kdz) to a global node number for the unknowns. Your function will look like nn=@(j,k,Nr) {some function of j, k, Nr}; (See Workshops 25-27 for examples of what this looks like). When setting up the matrix problem the function nn (or at least a similar process) can be used to reliably determine the row and column number of each entry in the Laplacian matrix and the row number in the BC matrix. Q3c Modify Lab06_Q3c.m to solve the 2D axisymmetric heat equation using the BTCS scheme (i.e. Implement your algorithm from Q3a) together with the function you have written in Q3b. Use the following values of the parameters • a = 0.05, • Nr = Nz =N (use 4 values 11, 21, 31 and 41), • time step Dt = 1 s • solve from t = 0 to t = 50 s for each value of N Use the commands ‘tic’ and ‘toc’ to determine how much time is spent setting the matrix problem up, (i.e. the time spent assembling the Laplacian matrix and BC vector) and separately, how much time is spent undertaking the time-stepping solution. (NOTE: The time taken for plotting should not be included as part of the time taken for time-stepping). Print out the time required for setting up each grid and solving each solution. Write a few sentences comparing the times and explain why you think you get the results you get. Using the subplot command, plot contours of your solution at time=50 including the boundary values for each (Nr, Nz) in a 2x2 grid. You should have one figure window displaying the contours as: MEC3456 Assignment Page 8 of 9 Note that you will have to reshape your 1D vector of unknowns into a 2D matrix and embed that into a slightly bigger 2D matrix so that you can add boundary conditions. Repeat this question for a time step of time step Dt = 50 s. Write a few sentences to the command window that compares the results for the 2 different timesteps. Are you surprised at your results for a single timestep of 50s compared to those for multiple steps with a 1 s timestep? QUESTION 4 [10 MARKS TOTAL] Q4a Suppose that the boundary condition is now a function of time, such: (, , ) = sin(2) 56.689 Eqn (6) Write a Matlab function that sets up the Laplacian matrix L and the boundary condition vector BC when the BCs vary in time. The only differences to Q3b will be the requirement to include time as an input parameter and the BC function will be as specified in Eqn 6. Your function header must be written as function [L,BC] = MatSet_T(Nr,Nz,rg,zg,tnow,alpha,BCf) where • tnow is the time at the current step, which is needed to evaluate the boundary condition function, • the rest of the input and output are the same as in Q3b. *Hint: You can modify from the function MatSet.m that you have written in Q3b, but you need to include both MatSet.m and MatSet_T.m as two separate MATLAB files. MEC3456 Assignment Page 9 of 9 Q4b Repeat Q3c with the new boundary condition as shown in Eqn (6). Instead of writing the time to set up and solve the system separately, print out the combined time for each different grid. Print to the command window a few sentences comparing the times in Q3c and Q4b, and explain what causes the difference. Compare the solutions between Q3 and Q4 - are there any differences? What are they and why? Write a few sentences to the command window to answer. THERE IS NO BONUS QUESTION FOR THIS LAB Poor Programming Practices [-10 Marks] (Includes, but is not limited to, poor coding style or insufficient comments or unlabeled figures, etc.) (END OF LAB)





























































































































































































































































































































































学霸联盟