Modelling and Simulation 2 (Dr M Jabbari): CFD CW 2a Due Thursday, April 30, 2020 – 17:00 CFD CW 2A – MATLAB Background Practical problems of heat transfer by conduction are often quite intricate and cannot be solved by ana-lytical methods. Their mathematical models may include non-linear differential equations with complexboundary conditions. In such cases, the only recourse is to obtain approximate solutions by employingdiscrete numerical techniques like the Finite Difference Method (FDM). A description of the physical problem is shown in Figure 1. A suitable model for this physical problemcould be based on Fourier’s law of heat transfer. Conduction through a thin layer of fluid (L << W and L << H) can be approximated in one dimension. The rate of heat transfer qx is proportional to negativetemperature gradient, qx = −kA∂T ∂t = −k (dydz) ∂T ∂t (1) H W L x y z i = 1 i = 2 i− 1 i i+ 1 i = N − 1 i = N x = 0 x = L ∆x −k dTdx ∣∣Left −k dTdx ∣∣Right Control volume boundary Figure 1: A 3D solid wall and its 1D FD discretisation along thickness. Now consider the energy balance for a small element “i” of fluid, shown in Figure 1, to be rate of heat conduction rate of heat conductioninto control volume out of control volume+ = +rate of heat generation rate of energy storageinside control volume inside control volume which can be expressed algebraically as qx + qgen = qx+∆x + qs (2) qx+∆x can then be expressed as a Taylor series expansion where, neglecting terms involving the secondderivative and higher, qx+∆x = qx + ∂qx ∂x dx Eq.1−−−→ qx+∆x = qx + ∂ ∂x ( −k∂T ∂x ) dxdydz (3) In the second half of Eq. 3 we have substituted in Eq. 1 for qx. The energy storage term may be expressedas Page 1 of 5 Modelling and Simulation 2 (Dr M Jabbari): CFD CW 2a Due Thursday, April 30, 2020 – 17:00 qs = ρcp ∂T ∂t dxdydz (4) where ρcp ∂T∂t is the time rate of change of the sensible (thermal) energy of the medium per unit volume.Substituting Eqs. 3 and 4 into Eq. 2 and neglecting any heat generation (like chemical reaction, qgen = 0),the energy balance equations becomes qx + 0 = qx + ∂ ∂x ( −k∂T ∂x ) dxdydz + ρcp ∂T ∂t dxdydz (5) Assuming a constant thermal conductivity, k, we can take it outside the differential operator. We can alsocancel the qx term and factor out the volume (dxdydz), giving us the final equation ρcp ∂T ∂t = k ∂2T ∂x2 ⇔ ∂T ∂t = k ρcp ∂2T ∂x2 ⇔ ∂T ∂t = α ∂2T ∂x2 (6) whereα = k/ρcp is the thermal diffusivity. In the presence of internal heat source/sink, Q˙, Eq. (6) becomes: ∂T ∂t = α ∂2T ∂x2 + Q˙ ρcp (7) Getting Started Eq. 7 is the model we are going to use for MATLAB programming. It is a two-dimensional problem, varyingover length (x) and time (t), c.f. Figure 2. With that in mind, we must construct a grid of points in spaceand a second grid of points in time. The spacing of these points are dictated by our chosen discrete timestep “∆t” and an analogous discrete spatial step “∆x”. The grid of points in both space and time can bedefined as follows tn = n∆t , n = 0, 1, ..., N xi = i∆x , i = 0, 1, ..., P (8) Nodal temperatures now depend on two indices, “i” and “n”, which correspond to the spatial and timepoints respectively: Tni ≡ T (xi, tn) (9) Figure 2: 2D distretisation domian. Page 2 of 5 Modelling and Simulation 2 (Dr M Jabbari): CFD CW 2a Due Thursday, April 30, 2020 – 17:00 Instructions In order to develop a FDM solver in Matlab, you will need to : • Develop finite difference approximations (from Taylor series) for the differential operators on theLHS and RHS of Eq. 7. • Substitute these expressions to produce a discrete finite difference equation in terms of the tem-perature at the grid points as illustrated in Figure 2. • Assume you are using a grid of P uniformly-spaced points in space (i = 1, ..., P ). • Rearrange this equation such that you have an expression for finding the temperature at the nexttime step in the simulation in terms of values at the previous time step. Hint 1: The LHS of the equation features a first derivative which can be discretised using either “first-orderbackward difference scheme” or “first-order forward difference scheme”. Hint 2: The derivative in the RHS features a second derivative which can be discretised using a “second-order central difference scheme”. You may need to combine Taylor series to generate the correct expres-sion. Hint 3: The source/sink term can be treated as a scalar. Page 3 of 5 Modelling and Simulation 2 (Dr M Jabbari): CFD CW 2a Due Thursday, April 30, 2020 – 17:00 Question 1 Question Weighting: 60% Task 1: In this part of the assignment, you will implement your equations in MATLAB, thus creatingyour very own 1D heat equation solver based on Eq. (7) and assuming that the source termis only applied in the centre-point (e.g. at i = 4 when P = 7). You will have to submit yourMATLAB code (and any accompanying functions) asm-file. Make sure you name your file in thisway: “StudentID CW2a.m”. Remember to avoid hard-coding values like the number of grid points P . It should be easy for a user of your code to change a setting in one place and that change isautomatically propagated throughout your code. Use the following information as global settings: 1 c l ea r a l l 2 c l c 3 4 %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− 5 % Input parameters 6 %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− 7 P = 7 ; % Numbers of nodes i n x−d i r e c t i o n 8 L = 1 ; % length [m] 9 k = 40 ; % Thermal conduc t i v i t y [W/mK] 10 rho = 1000; % Dens i ty [ kg /m^3 ] 11 cp = 100 ; % Spec i f i c heat [ J /kgK ] 12 t f = 100 ; % Tota l t ime [ s ] 13 dt = 1 ; % dt va lue [ s ] 14 T _ l e f t = 50. ; % T at l e f t boundary [ C ] 15 T_ r i gh t = 100. ; % T at r i g h t boundary [ C ] 16 Q = 500000. ; % Source term [W/m^3 ] 17 %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− Initialise your grid points with and initial temperature of 20 ◦C i.e. for i = 1 : P ⇒ T (xi, t0) = 20 ◦C (10) For boundary conditions at both ends of the spatial domain use: T (x0, tn) = Tleft = 50 ◦C T (xP , tn) = Tright = 100 ◦C (11) Hint 1: Use a while loop to march your simulation through time, solving your discrete equation foreach time step. Hint 2: You may want to use the code layout given below: 1 %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− 2 %% Time loop 3 whi le ( time <= t f ) 4 %increment time 5 6 %solve fo r T 7 8 %Ex t r ac t the temperature and put i n the T g r i d 9 end 10 11 %% p l o t t i n g 12 f i g u r e 13 %−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− Question 1 continued on next page. . . Page 4 of 5 Modelling and Simulation 2 (Dr M Jabbari): CFD CW 2a Due Thursday, April 30, 2020 – 17:00 Task 2: Use MATLAB to plot the temperature along x-direction (T − x diagram) at the end of 100seconds — (tagged: Fig. 1) Task 3: Use MATLAB to plot the temperature at each time step at grid points “2”, “4” and “6”, T (x2, tn), T (x4, tn) and T (x6, tn), over time (T − t diagram) — (all in one graph and tagged: Fig.2) Question 2 Question Weighting: 40% Use the program you developed in Question 1 to run two different tasks: Task 1: Run the simulations with different grid resolution, P = 9, 11, 29, and compare the “T − x”at t = 100s (all in one graph and tagged: Fig. 3) and “T − t” diagram for the centre grid point(all in one graph and tagged: Fig. 4). Discuss your results briefly in terms of numerical modellingaccuracy. Specifically, how does changing the resolution of the grid change the solution? Task 2: Run the simulations with grid resolution of P = 29, but different final time of tf = 100, 400, 700, 1000s and compare “T − x” diagrams as well as “T − t” diagram (all in one graphand tagged: Fig. 5) for the centre grid point (all in one graph and tagged: Fig. 6). Discuss yourresults briefly in terms of steady/unsteady phenomena. Is this problem unsteady or steady? Whathappens when you run your simulation for longer? Assessment You will need to submit yourm-file and a single PDF document in the blackboard that includes Figs. 1–6 and all necessary discussions. This assessment is worth 10% of the whole unit. The marking scheme for this coursework is available in the blackboard (marking scheme CW2.pdf). Page 5 of 5
