SOLUTION: dimensional convection diffusion | My Assignment Tutor

Computational Fluid Dynamics Chapter 031Chapter 3One-dimensional convection diffusion problemsWhen we study moving of heat or a substance in fluid flow, apart from diffusion, convectionalso makes contribution, and it is very common that convection makes more contribution thandiffusion. Convection is heat transfer by mass motion of a fluid such as air or water when the heatedfluid is caused to move away from the source of heat, carrying energy with it. In this chapter, wediscuss the FDM solution of the one-dimensional convection-diffusion equation.2.1 One-dimensional (1D) convection-diffusion equationWe consider a flow in a tube as shown in Figure 1. The fluid flows in the tube with a velocity u.If the temperature of the fluid is not evenly distributed in the tube, there should be diffusion of theheat in the tube. The heat also flows with the fluid flow. The temperature of the fluid in the tube isinitially 20°C and then the 80°C water starts to flow into the tube from the left end. In this example,the initial and boundary conditions are:Initial condition:T(x)  20C , throughout the tube (1)Boundary condition:80 CT x0   ,  0 Tx xL(2)whereL=the length of the tubeu = the fluid velocity in the pipeFigure 1 One-dimensional convection and conduction of heat in a tubeThe heat will be transferred towards the right end of the tube due to the convection (going withthe flow) and diffusion (going from a high temperature area to a low temperature area). We choose avery small segment of the tube Δx. The heat flow rate (energy passing unit area in unit time) at theleft and right sides of Δx are q1 and q2, respectively. The heat flow rate is related to the velocity andthe temperature asΔxq1 q2Flow velocity in pipe is uxx=0 x=LComputational Fluid Dynamics Chapter 032Txq ucpT K    (3)The change rate of the temperature of Δx is related to q as (see Chapter 2)qxt cTp   1(4)Substituting Equation (4) into Equation (3), we get 22TxxuTTt    (5)wherecpK  (6)If the velocity u is constant, then Equation (5) becomes22TxTxuTt     (7)Equation (7) is the convection-diffusion equation for solving the problem defined in Figure 1.The termTxu  is the convection term and22Tx  is the diffusion term.2.2 Solution of the one-dimensional (1D) convection-diffusion equation using the ForwardTime Central-Space Scheme (FTCS)Now we solve Equation (5) using FDM. We divide the tube into many segments. The length ofeach segment is Δx as shown in Figure 2.Figure 2 Computational grid for the convection diffusion problemIf we replace the derivatives in Equation (7) using the following FDM approximations (they arederived in the previous chapter)x=L xx2=Δxxx1=0x3=2 Δxxi=(i-1) Δxxm+1=mΔxuComputational Fluid Dynamics Chapter 033 21 12xT T Tninini n i 2x2T   (8)Tx Tnin1ni1i2 Tx  (9)tT T Tt ninin i 1(10)we can get2 1Tu1tT T 1 1 12T T Tn  2 xTxinininininini   (11)Then, the FDM formula for calculating the temperature next step becomesnininini TxttxuTxtTxttxuT2 1 2 2 1121 22                 (12)If we definetxua  and2xtb (13)Equation (12) becomesTin1    a  b Tin 1  1 2bTin   a  bTin1 (14)Equation (14) is the Forward-Time Central-Space Scheme (FTCS).Appendix I is an example program for solving the convection equation. C in the program standsfor temperature. The problem is the heat and fluid flow in a tube as discussed in the beginning ofthis chapter. Figure 3 is an example of the solution of the convection-diffusion equation (Equation(7)).The parameters in the program in Appendix I are: Velocity u=0.1; Diffusion coefficient α=0.001.The computational grid is: Length of the computational domain is 5 Grid number = 100 Grid size Δx = 0.05 Time step Δt=0.1The initial and boundary conditions of this problem are: Initial condition: the temperature is zero along the whole tube Boundary condition: (a) the temperature at the left end is 1; (b) at the right end, thetemperature gradient in the x-direction is zero.Computational Fluid Dynamics Chapter 034Figure 3 is the temperature distribution in the tube after 200 time step (at t  20  200 0.1). Itcan be seen that the numerical results follow the trend of the analytical solution very well, althoughthere is some error in the numerical results at x=1.5 – 2.Figure 3 Distribution of the temperature (here C represents the temperature) along the x-directionobtained by solving the convection diffusion equation using the program in Appendix I (FTCSmethod).2.3 Solution of one-dimensional (1D) convection-diffusion equation using the implicit methodForward-Time Central-Space Scheme (FTCS) is an simple explicit method, the first and secondorder spatial derivatives of the temperature T in Eq. (7) are calculated using the values at the timent , which has been known.If the spatial derivatives are calculated using the values at the time (n 1)t , which are notknown, the scheme is implicit method. The FDM discretization of the equation is:211 1 112T T Tn 1T TTu (15)i 11112 xTxtnininininini By moving all the values at time (n 1)t to the left hand side of the equation, Eq. (15)becomesnininini T TxttxuTxtTxttxu             12 11212 121 22  (16)Eq. (16) can be rewritten asnini ini in aiTi1 bT c T  T 1 1 11 (17)where0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 500.20.40.60.811.21.4xCnumerical,dt=0.1analyticalComputational Fluid Dynamics Chapter 035ai   2uxt  x2t ,bi 1 2x2t ,ci  2uxt  x2t (18)Equation (17) is the implicit FDM formula for the convection-diffusion equation. At every FDMgrid, there is an Equation (17). By putting all the equations together, we have a matrix equation:       Nn Nn n n nn nNnNN NN N Nn n ndd d d d dT T T T T Ta ba b ca b ca b ba b cb c11 2 31 1 111312111 1 13 3 32 2 21 1  (19)Because each row of the coefficient matrix in Equation (19) has only three nonzero elements andall the nonzero elements aligned diagonally in the matrix, the coefficient matrix in Eq. (19) is calledtridiagonal matrix. In the computer program, we may use a one-dimensional array A, B, C and D tosave the coefficients. The tridiagnal equation can be solved analytically.In Appendix II, there is an example Matlab program for solving the 1D convection-diffusionequation using the implicit method. The function “solve_diagonal_eq” in Appendix II is the Matlabfunction for solving the tridiagnal equation.Figure 4 is an example of the numerical solution of the 1D convection-diffusion equation usingthe implicit FDM. The calculation parameters are the same as those used for Figure 3. It can be seenthe curve generated by the implicit method is smoother than that by the FTCS method.Figure 4 Solution of the convection-diffusion equation using the implicit FDM.0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5-0.200.20.40.60.811.2xTnumerical,dt=0.1analyticalComputational Fluid Dynamics Chapter 036Appendix IA MATLAB program for solving 1D convection-diffusion equation using the FTCS scheme,including the main program and a function for calculating the analytical solution.%A MATLAB program for solving 1D convection-diffusion equation%using Forward-time Central-space scheme (FTCS)% u — the velocity% D — the diffusion coefficient% C — the density to be calculated, or the temperature% t — time% x — x-coordinate array for all grids%% dx– spatial time step% dt– temperal time step% N — total time steps to be calculated; final time t = N*dt% L — length of domain% m — grid number is m+1;% Initial condition C=0 at whole domain except x=0 where C=C0% boundary condition C=C0 at x=0% C1– saves C at previous time step% ya– analytical solution of C% p03_1 — function for calculating analytical solutionclear;u = 0.1;C0 = 1;D = 0.001;L = 5;m = 100;dx = L/m;x = (0:m)*dx;t = 0;dt = 0.1;N = 200;a = u*dt/(2*dx);b = D*dt/(dx*dx);C(1:m+1) = 0;C(1)= C0;for n=1:Nt = t+dt;C1 = C;for i=2:mC(i) = (a+b)*C1(i-1)+(1-2*b)*C1(i)+(-a+b)*C1(i+1);endC(m+1) = C(m);d =1; % USELESSendplot(x,C,’o’);xlabel(‘x’); %honrizontal axis’s label is ‘x’ylabel(‘C’); %vertical axis’s label is ‘T’Computational Fluid Dynamics Chapter 037 hold on;ya = p03_1(x, t, C0, u, D,L);plot (x,ya,’r’);h = legend(‘numerical,dt=0.1′,’analytical’);%MATLAB program for analytical solution of 1D convection-diffusionequation% u — the velocity% D — the diffusive coefficient% C — the density to be calculated% t — time% x — x-coordinate array for all grids% Initial condition C=0 at whole domain, except x=0 where C=C0% boundary condition C=C0 at x=0% C1– saves C at previous time stepfunction C = p03_1(x, t, C0, u, D,L)A1 = 0.5*erfc((x-u*t)/(2*sqrt(D*t)));A2 = 0.5*exp(u*x/D) .*erfc((x+u*t)/(2*sqrt(D*t)));A3 = 0.5*(2+u*(2*L-x)/D+u^2*t/D);A4 = exp(u*L/D)*erfc(((2*L-x)+u*t)/(2*sqrt(D*t)));A5 = (u^2*t/(acos(-1)*D))^2.*exp(u*L/D-(2*L-x+u*t).^2/(4*D*t));A = A1+A2+A3.*A4+A5;C = C0*A;end Appendix IIA Matlab program for solving the convection-diffuson equations using implicit method%MATLAB program for solving 1D convection-diffusion equation%using Forward-time Central-space scheme (FTCS)% u — the velocity% Diff — the diffusive coefficient% T — the density to be calculated% t — time% x — x-coordinate array for all grids%% dx– spatial time step% dt– temperal time step% N — total time steps to be calculated; final time t = N*dt% L — length of domain% m — grid number is m+1;% Initial condition T=0 at whole domain except x=0 where T=C0% boundary condition T=C0 at x=0% T1– saves T at previous time step% ya– analytical solution of T% analytical — function for calculating analytical solution% solve_diagonal_eq — function for solving the tridiagonal matrixeqComputational Fluid Dynamics Chapter 038clear;u = 0.1;C0 = 1;Diff = 0.001;L = 5;m = 100;dx = L/m;x = (0:m)*dx;t = 0;dt = 0.1;N = 200;T(1:m+1) = 0;A(1:m+1) = 0;B(1:m+1) = 0;C(1:m+1) = 0;D(1:m+1) = 0;T(1)= C0;for n=1:Nt = t+dt;T1 = T;A(1) = 0;B(1) = 1;C(1) = 0;D(1) = C0;A(m+1) = 1;B(m+1) = 1;C(m+1) = 0;D(m+1) = 0;for i=2:mA(i) = -u*dt/(2*dx) – Diff*dt/(dx*dx);B(i) = 1+2*Diff*dt/(dx*dx);C(i) = u*dt/(2*dx) – Diff*dt/(dx*dx);D(i) = T(i);endM=m+1;T=solve_diagonal_eq(M,A,B,C,D);endplot(x,T,’o’);xlabel(‘x’); %honrizontal axis’s label is ‘x’ylabel(‘T’); %vertical axis’s label is ‘T’hold;ya = analytical(x, t, C0, u, Diff,L);plot (x,ya,’r’);Computational Fluid Dynamics Chapter 039 h = legend(‘numerical,dt=0.1′,’analytical’);% function for solving a tridiagonal matrix equation% M is the dimension of the arrays% A, B, C are the triadiagonal matrix coefficients% D is the vector at the right of the equation% T is the solution vectorfunction T=solve_diagonal_eq(M,A,B,C,D)C(1) = -C(1)/B(1);D(1) = D(1)/B(1);TERM = 0;for I=2:MTERM = B(I)+A(I)*C(I-1);TERM = 1.0/TERM;C(I) = -C(I)*TERM;D(I) = (D(I)-A(I)*D(I-1))*TERM;endfor I=M-1:-1:1D(I) = C(I)*D(I+1)+D(I);endT=D;end