SOLUTION: Computational Fluid Dynamics | My Assignment Tutor

Computational Fluid Dynamics Chapter 21Chapter 2 One-dimensional diffusion problems2.1 One-dimensional (1D) diffusion equationDiffusion problems are very common in engineering. Diffusion is the movement of a substancefrom a region of higher concentration to a region of lower concentration. For example, the heattransfer from high temperature area to low temperature area is due to conduction. In numericalstudies, conduction is solved by the diffusion equation. Another example is the pollutant in theatmosphere moving from the high density to low density areas through diffusion.Consider a metal rod shown in Figure 1, initially the temperature of the whole rod is 0°C. Nowthe temperature at the two ends of the rod suddenly becomes 100°C and remains to be 100°C for therest of the time. We need to find out how the temperature is changing along the whole rod. Weassume the heat only transfer along the axial direction of the rod and there is not any heat exiting orentering the rod through the rod surface.The problem defined in Figure 1 is a typical one-dimensional (1D) heat convection problem.Heat will be transferring from the two ends towards the center of the rod due to the temperaturedifference. When we solve the unsteady problem, we need both initial and boundary conditions. Theinitial and boundary conditions of the problem areT x Ct ( ) 0 0 (1-1)T 100C, at x  0 and x 1m (1-2)Figure 1 One dimensional heat diffusion problem. Initially the temperature of the whole rod is 0°C.The temperatures at the two ends become 100°C instantly and remain to be 100°C.To facilitate the analysis, the following parameters are defined: q = heat flux, the energy through a unit area of the cross-section of the rod in unit time cp=specific heat, the heat (energy) required to increase the temperature of unit mass of the rodfor one degree ρ=density of the rod material A=cross sectional area of the rod T=temperature ΔT=change of the temperature Δt=time step, small period of timeBased on the heat transfer theory, the heat flux (the energy therough a unit area of the crosssection) is proportional to the temperature gradient, i.e.Txq K   (2) Metal rodThermal insulation surfacexx=0x=1 mT = 100°CT = 100°C Computational Fluid Dynamics Chapter 22If we choose a small segment of the rod with a length of Δx, the heat transfer at the two ends ofthe small segment Δx is q1and q2, respectively, as shown in Figure 2. If the temperature of the rod Δxincreases by ΔT, the required heat isFigure 2 Heat transfer in a very small segment of the rod with a length of Δx c m T c A x TEnergy required  p   p  (3)The energy through the left and right ends of Δx are q1 and q2, respectively. If are q1 and q2 arenot the same, the net energy received by Δx isenergy gain  tA(q1  q2)(4) Based on Equation (4), if q1>q2, Δx receives heat and the temperature will increase and vice versa.Based on the energy balance equation (conservation of energy), the quantities in Equations (2) and (3)should be the same as each other, i.e. c pAxT  tA(q1  q2)Equation (5) can be rewritten as the change rate of the temperature:(5) xq qc1 2 1t  Tp  (6)If the length of the segment Δx is infinitesimal, Equation (6) becomes a partial differential equation(PDE)qxt cTp   1(7)Utilyzing Equation (2), i.e.Txq K   , Equation (7) becomes22TxKcTtp  (8)Equation (8) is the one-dimensional diffusion equation. If we definec pK  (9)xq1 q2Flow velocity in pipe is uxComputational Fluid Dynamics Chapter 23Equation (8) becomes22TxTt   (10)2.2 Solution of 1D diffusion equation using Finite Difference Method (FDM)Now we solve the 1D diffusion equation (10) numerically using Finite Difference Method(FDM). In FDM, the rod between x=0 and x=1 m is divided into many small segments as shown inFigure 3 (a). The length of each segment is Δx. The grid system shown in Figure 3 (a) is called thecomputational grid. Each location on the grid is called a node.When we solve the diffusion equation using FDM, we know the temperature distribution alongthe rod at t=0, i.e. the initial condition. Based on the initial condition, we can use the diffusionequation (10) to calculate the temperature along the rod after a very small time interval (time step) Δt.Then based on the temperature at t=Δt, we calculate the temperature at t=2Δt, 3Δt, …The temperature of the ith node (xi) at time t=nΔt is defined as Tin as shown in Figure 3 (b).(a)(b)Figure 3 (a) Spatial gride; (b) Temporal and spatial gridx=1 m xx2xxx1x3xxi(i-1)xXm+1mxxt0=0 1 x2 xi xmt1=ttn=(n)tnTiComputational Fluid Dynamics Chapter 24In FDM, the partial differential operators are approximated by the FDM formulae. Take a node i asan example as shown in Figure 4, the first order spatial derivativeTx is approximated byTx Tnin1ni1i2TTnininxixiT Tninin Tx  (11) Tx  1 (12) Tx  1 (13)Equation (11) is a central differential method and has second order accuracy. Equation (12) and(13) only has first order accuracy, less accurate than Equation (11).The first order spatial derivatives of T at i-1/2 and i+1/2 are TxTnini1ni Tx  12and xT Tnini1n i Tx  12(14)Base on Equation (14), the FDM formula for the second order derivative of T is21 11 1121222 2xT T TxTxTxT Tx Tx Tx Txnininininininin in i       (13)or 2xi2x2Tnini 1 12T T Tnin    (14)Equation (14) is the FDM formula for the second order derivative of T with respective to x.Figure 4 Definition of a node point inTi1nTinTi1ii-1 i+1xi-1/2 i+1/2xComputational Fluid Dynamics Chapter 25On the ith node, if we replacen i Tt by the first order FDM formulatT T nini 1, and replace thesecond order spatial derivativen i 2T2x by 1 2 2 1xT T T ninini   , the diffusion Equation (10) becomes21 11 2xT T TtT T ninininini    (15).If we know the values at the nth time step (at t=nΔt), based on Equation (15), the temperature atthe (n+1)th time step can be calculated by TinsT1  sT sTin1 1 in12i n (16)where2xa ts  (17)Equation (16) is the first order in time, second order in space, explicit FDM formula for the 1Ddiffusion equation. Table 1 is a MATLAB program code for solving the 1D diffusion equation (10).Figure 5 is the picture generated from the program. The boundary condition and the initial conditionof the problem are: Initial condition: temperature is 0°C for the whole the rod at t=0 Boundary condition: the temperature at the two ends of the rod is always 100°CIt can be seen that at t=1 s, the temperatures at the two ends are 100°C. Because the heat transferfrom the two ends to the center of the rod, the temperature at the center has increased from 0°C toabout 44°C. At t=5 s, the temperature at the center has been above 90°C. The temperature of thewhole rod will be close to 100°C after a long period of time.Table 1 A Matlab program for solving 1D diffusion Equation 10.clear; %clear all memorym = 16; %grid number in x is m+1dx = 1/m; %delta xx = (0:m)*dx; %coordinate at all grid pointsM = m+1; %Grid No. from 0 to ma = 0.1d0; %coefficient αt = 0; %present timedt = 0.01; %delta tN = 50; %solution is at time = N*dtT(1:m+1) = 0; %initial valueT(1) = 100; %initial at left endT(M) = 100; %initial at the right ends = a*dt/(dx^2);for k=1:N; %for each time stept = t+dt; %present time is tT0 = T; %T0 saves temperature at previous time stepComputational Fluid Dynamics Chapter 26T(1) = 100; %T at left boundary (end)T(M) = 100; %T at right boundary (end)for i=2:M-1 %to calculate temperature along the rodT(i) = T0(i)+a*dt/dx^2*(T0(i-1)-2*T0(i)+T0(i+1)); %for ecah grid points, calculate Tendendplot (x,T,’x’); %to plot out the x-T curve with symbol ‘o’xlabel(‘x’); %honrizontal axis’s label is ‘x’ylabel(‘T’); %vertical axis’s label is ‘T’hold on; %next plot will be in the same diagramxa = (0:200)/200; %xa isTa = p02_1(xa,a,t,20000); %Ta is analytical solution of Tplot(xa(1:201),Ta(1:201),’-r’); % plot out the Ta in the same diagram used for Th = legend(‘numerical’,’exact’);Table 2 A Matlab program for solving 1D diffusion Equation 10.MATLAB code for the analytical (exact) solution, the function p02_1 will be called by theprogram listed in Table 1File name: P02_1.mfunction y = p02_1(x, a, t, M) %analytical solutiony = 100;pi = acos(-1);for m=1:Mm1 = 2*m-1;y = y-400/(m1*pi)*sin(m1*pi*x)*exp(-a*m1^2*pi^2*t);end(a) t=1 s (b) t= 5 sFigure 5 FDM solution of the 1D diffusion equation using the program listed in Table 1.0 0.2 0.4 0.6 0.8 15060708090100xTnumericalexit Exact0 0.2 0.4 0.6 0.8 19999.299.499.699.8100xTdx=0.125,dt=0.05dx=0.0625,dt=0.01exactComputational Fluid Dynamics Chapter 27In this section, we introduced an explicit FDM formula for solving the 1D diffusion equation. Ithas been seen that the accuracy of the solution depends on the grid size Δx and time step Δt. Whensolving the initial-value problem, the numerical scheme should be stable. For a particular grid size,the solution is not stable if the time step is too big, the scheme will be unstable as shown in Figure 6,where the meaningless values are obtained because the time step is too big. The stable condition for anumerical scheme will be discussed in the future chapters.Figure 6 If the time step is too big, the solution is unstable-8.E+33-6.E+33-4.E+33-2.E+330.E+002.E+334.E+336.E+338.E+330 0.2 0.4 0.6 0.8 1Temperature TxΔx=0.05 m, Δt=0.02 s