Chapter 32 According to Fick’s First Law, the Diffusive Flux is

CHAPTER 32 32.1 First equation 6.075c 0 − 3.2 c1 = 262.5 Middle equations (i = 1 to 8) −2 .1c i −1 + 3.45c i − 11 . c i+1 = 0 Last equation −3.2 c 8 + 3.45c 9 = 0 The solution is 32.2 Element equation: (See solution for Prob. 31.4 for derivation of element equation.)  a 11   a 21 a 12   c1   b1    =   a 22   c 2  b2  a 11 = 2 1 0.2 − + ( 2 .5) = 0.55 2.5 2 2 a 12 = −2 1 + = −0.3 2.5 2 a 21 = −2 1 − = −1.3 2.5 2 a 22 = 2 1 0.2 + + ( 2.5) = 155 . 2 .5 2 2 b1 = −2 dc (x ) dx 1 b2 = 2 Assembly:  0.55 −0.3  c 0     −1.3 2 .1 −0.3 c    1    c 2  = −1.3 2 .1 −0.3   −1.3 2 .1 −0.3 c 3      −1.3 155 .  c 4  Boundary conditions: Inlet: dc ( 0) dx dc Uc 0 − Uc in ( 0) = dx D Uc in = Uc 0 − D Substitute into first equation  dc  − dx ( x 1 )    0   0     0  dc   (x 2 )   dx  dc (x ) dx 2 0.55c 0 − 0.3c1 = 100 − c 0 1.55c 0 − 0.3c1 = 100 Outlet: dc (10 ) = 0 dx Solution: c0 = 74.4 c1 = 51.08 c2 = 35.15 c3 = 24.72 c2 = 20.74 32.3 According to Fick’s first law, the diffusive flux is J( x) = −D dc ( x) dx where J(x) = flux at position x. If c has units of g/m3, D has units of m2/d and x is measured in m, flux has units of g/m2/d. In addition, there will be an advective flux which can be calculated as J ( x ) = Uc ( x ) Finite divided differences can be used to approximate the derivatives. For example, for the point at the beginning of the tank, a forward difference can be used to calculate dc dx ( 0) ≅ 52.47 − 76.44 2.5 = −0.9588 g / m3 m Thus, the flux at the head of the tank is J ( x ) = −2 ( −0.9588 ) + 1( 76.44 ) = 19 .176 + 76.44 = 95.616 g / m3 m The remasinder of the values can be calculated in a similar fashion using centered (middle nodes) and backward differences (the end node): 32.4 Segmentation scheme: 1 ,2 2 ,2 3 ,2 4 ,2 5 ,2 6 ,2 1 ,1 2 ,1 3 ,1 4 ,1 5 ,1 6 ,1 1 ,0 2 ,0 3 ,0 40 40 40 100 d c /d y = 0 100 d c/d x = 0 d c /d y = 0 100 Nodes 1,1 through 5,1 0.4 c i +1, j − 2 ci , j + c i −1, j 52 + 0.4 c i , j +1 − 2 ci , j + c i , j +1 52 − 0.2 ci , j Collecting terms gives 0.264 ci , j − 0.016 ci +1, j − 0.016ci −1, j − 0.016 ci , j +1 − 0.016c i , j +1 = 0 Node 6,1 would be modified to reflect the no flow condition in x and the Dirichlet condition at 6,0: 0.264 c 6,1 − 0.032 c 5,1 − 0.016c 6, 2 − 0.016 (100 ) = 0 The nodes along the upper edge (1,2 through 5,2) would be generally written to reflect the no-flow condition in y as 0.264 c i , j − 0.016c i +1, j − 0.016c i −1, j − 0.032 ci , j +1 = 0 The node at the upper right edge (6,2) would be generally written to reflect the no-flow condition in x and y as 0.264 c 6, 2 − 0.032 c 5,2 − 0.032 c 6,1 = 0 Finally, the nodes along the lower edge (1,0 through 3,0) would be generally written to reflect the no-flow condition in y as 0.264 c i , j − 0.016c i +1, j − 0.016c i −1, j − 0.032 ci , j +1 = 0 These equations can be solved for 32.5 For simplicity, we will use a very coarse grid to solve this problem. Thus, we place nodes as in the following diagram. x v10o oil v8o v6o v6w v4w v2w 0 water v0w A simple explicit solution can be developed by substituting finite-differences for the second derivative terms in the motion equations. This is done for the three non-boundary nodes, dv 2 w v − 2v 2 w + v 4 w = µw 0 w dt ∆x 2 dv 4 w v − 2v 4 w + v 6 w = µ w 2w dt ∆x 2 dv 8o dt = µo v 6o − 2 v8o + v 10o ∆x 2 These three equations have 7 unknowns (v0w, v2w, v4w, v6w, v6o, v8o, v10o). The boundary conditions at the plates effectively specify v0w = 0 and v10o = 7. The former is called a “no slip” condition because it specifies that the velocity at the lower plate is zero. The relationships at the oil-water interface can be used to used to eliminate two of the remaining unknowns. The first condition states that v6o = v6w (i) The second can be rearrange to yield v 6w = µ o v 8o + µ w v 4 w (ii) µo + µw These, along with the wall boundary conditions can be substituted into the differential equations dv 2 w − 2v 2 w + v 4 w = µw dt ∆x 2 dv 4 w dt v 2 w − 2v 4 w + = µw µ oν 8o + µ w ν 4 w ∆x µo + µw 2 dv 8o dt = µo µ oν 8o + µ wν 4 w − 2 v 8o + 7 µo + µw ∆x 2 These equations can now be integrated to determine the velocities as a function of time. Equations (i) and (ii) can be used to determine v6o and v6w. The results are plotted below: 8 t = 1.5 s t = 1.0 s 4 t = 0.5 s 0 0 2 4 6 8 10 32.6 Using a similar approach to Sec. 32.2, the following nodal equation can be developed for node 11: 4 u 11 − 1.21954 u12 − 1.21954 u 10 − 0.78049 u12 − 0.78049 u 01 = 0.357866 Similar equations can be written for the other nodes and the resulting equations solved for A graphical comparison of the results from Sec. 32.2 can be made with these results by developing a plot along the y dimension in the middle of the plate: 0 0.5 1 1.5 2 0 -0.3 -0.6 These results can then be used as input to the right-hand side of Eq. 32.14 and the resulting simultaneous equations solved for Again the comparison is good 0.12 0.08 0.04 0 0 0.5 1 1.5 2 32.7 Grid scheme 0, 2 1, 2 2, 2 3, 2 0, 1 1, 1 2, 1 0, 0 1, 0 2, 0 3, 1 All nodes in the above scheme can be modeled with the following general difference equation hi +1, j − 2 hi , j + hi −1, j ∆x 2 + hi , j +1 − 2 hi , j + hi , j −1 ∆y 2 =0 Node 0,0: h1, 0 − 2 h0 ,0 + h−1, 0 ∆x 2 + h0,1 − 2 h0, 0 + h0, −1 ∆y 2 =0 The external nodes can be approximated with finite differences dh h 0,1 − h0, −1 = dy 2∆y h0, −1 = h0,1 − 2∆y dh = h0,1 dy dh h1, 0 − h−1,0 = dx 2∆x h−1,0 = h1,0 − 2∆x dh = h1,0 − 2(1)(1) = h1,0 − 2 dx which can be substituted into the difference equation to give 2h1,0 − 2h0,0 − 2 2 4h0, 0 + 2h0,1 − 2h0 ,0 ∆x ∆y 2 − 2h1, 0 − 2h0 ,1 = −2 =0 Node 1,0: 4 h1,0 − 2 h1,1 − h0 ,0 − h2 ,0 = 0 Node 2,0: 4 h 2, 0 − 2 h1,0 − 2h 2 ,1 = 0 Node 0,1: 4h0,1 − 2h1,1 − h0,0 − h0,2 = −2 Node 1,1: 4 h1,1 − h1, 0 − h 0,1 − h1,2 − h2 ,1 = 0 Node 2,1: 4 h 2,1 − h1,1 − h2 ,2 − h3,1 − h2 ,0 = 0 Node 0,2: 4h0 ,2 − 2h0,1 − 2h1, 2 = −2 Node 1,2: 4h1,2 − h0,2 − h2 ,2 − 2 h1,1 = 0 Node 2,2: 4h2 ,2 − h1,2 − 2 h2,1 = 20 The equations can be solved simultaneously for More refined results can be obtained by using a finer grid spacing. 32.8 The fluxes can be determined using finite divided differences as 32.9 Because of the equi-spaced grid, the domain can be modeled with simple Laplacians. The resulting solution is 32.10 A convenient segmentation scheme can be developed as Simple Laplacians reflecting the boundary conditions can be developed and solved for 32.11 The system to be solved is  2.7  −2     −2   x 1   0   x   0 2.75 −0.75   2  =   −0.75 2.25 −1.5  x 3   0  − 15 . 1.5   x 4  2  which can be solved for x1 = 2.857, x2 = 3.857, x3 = 6.5238, and x4 = 7.857. 32.12 The system to be solved is −0.4 1.8  0.6  −0.4      −1.4 −1.4 2.1 −0.7 −0.7 1.6 −0.9   x 1   0   x   0   2      x 3  =  0  −0.9  x 4   0     0.9   x 5   1 which can be solved for x1 = 5, x2 = 7.5, x3 = 8.214286, x4 = 9.64286, and x5 = 10.75397. 32.13 Substituting the Crank-Nicolson finite difference analogues to the derivatives ∂ 2u ∂x 2 ∂u = ∂t 1  u i+1, n+1 − u i ,n +1 + u i−1, n+1 u i+1,n − u i ,n + u i−1,n  +   2 ∆x 2 ∆x 2  u i, n +1 − u i, n = ∆t into the governing equations gives the following finite difference equations:    ∆x 2 u + [ 1 ] u = − u + 2 − 2   i +1, n +1 i −1, n ∆t  i, n+1 ∆t  [1]u i −1, n+1 + − 2 − 2 ∆x   2  u i, n − ui +1, n  0≤x≤ 1 2   ∆x 2  1 ≤ x ≤1 u i, n+1 + [ r ]u i+1,n +1 = − ru i −1, n +  2r − 2 u i, n − ru i +1, n ∆t  ∆t  2   Substitute for the end point boundary conditions to get the end point finite difference equations. Substitute the first order Crank Nicolson analogues to the derivatives [ r ]u i−1, n+1 +  − 2r − 2 ∆ x 2 ∂u 1  u i +1, n+1 − ui −1, n+1 u i +1,n − u i −1, n  into the midpoint boundary condition and get =  +  ∂r 2  2∆r 2 ∆r  u a L +1,n +1 + u a L +1,n + r (u b L +1, n+1 + u b L +1, n ) = u L −1,n +1 + u L =1,n + r (u L +1, n+1 + u L +1, n ) where u a and u b are fictitious points located in the opposite side of the midpoint from their half. Write out the two finite difference equations from above for the point i = L (the midpoint) then combine these two equations with the midpoint boundary condition to obtain the midpoint finite difference equation:    ∆x 2  u L , n+1 + [ 2]u L +1,n +1 = −2u L −1, n + 2 (1 + r ) − 4 u L , n − (1 + r )u L +1,n ∆t  ∆t   [ 2]u L−1, n+1 + − 2 (1 + r ) − 4 ∆x  2 %PDE Parabolic Problem - Transient Heat conduction in a composite rod % u[xx]=u[t] 0< x< 0.5 % r(u[xx])=u[t] 0.5< x< 1 % BC u(0,t)=1 u(1,t)=1 % u[x]=r(u[x]) x=0.5 % IC u(x,0)=0 0< x < 1 % i=spatial index, from 1 to imax % R = no. of x points (R=21 for 20 dx spaces) % n=time index from 1 to N % N = no. of time steps, % Crank-Nicolson Formulation R=41; %(imax must be odd for point L to be correct) N=69; % last time step = nmax+1 L=(R-1)/2+1; % L = midpoint of point no. (for R=41, L=21) % Constants r=0.01; dx=1/(R-1); dx2=dx*dx; dt=dx2; % Setting dt to dx2 for good stabilility and results % Independent space variable x=0:dx:1; % Sizing matrices u=zeros(R,N+1); t=zeros(1,N+1); a=zeros(1,R); b=zeros(1,R); c=zeros(1,R); d=zeros(1,R); ba=zeros(1,R); ga=zeros(1,R); up=zeros(1,R); % Boundary Conditions at t=0 u(1,1)=1; u(R,1)=1; % % Time step loop n=1 represents 0 time, next time = n+1 t(1)=0; for n=1:N t(n+1)=t(n)+dt; % Boundary conditions & Constants u(1,n+1)=1; u(R,n+1)=1; dx2dt=dx2/dt; % coefficients b(2)=-2-2*dx2dt; c(2)=1; d(2)=(2-2*dx2dt)*u(2,n)-u(3,n)-2; for i=3:L-1 a(i)=1; b(i)=-2-2*dx2dt; c(i)=1; d(i)=-u(i-1,n)+(2-2*dx2dt)*u(i,n)-u(i+1,n); end a(L)=2; b(L)=-2*(1+r)-4*dx2dt; c(L)=2*r; d(L)=-2*u(L-1,n)+(2*(1+r)-4*dx2dt)*u(L,n)-2*r*u(L+1,n); for i=L+1:R-2 a(i)=r; b(i)=-2*r-2*dx2dt; c(i)=r; d(i)=-r*u(i-1,n)+(2*r-2*dx2dt)*u(i,n)-r*u(i+1,n); end a(R-1)=r; b(R-1)=-2*r-2*dx2dt; d(R-1)=-r*u(R-2,n)+(2*r-2*dx2dt)*u(R-1,n)-2*r; % Solution by Thomas Algorithm ba(2)=b(2); ga(2)=d(2)/b(2); for i=3:R-1 ba(i)=b(i)-a(i)*c(i-1)/ba(i-1); ga(i)=(d(i)-a(i)*ga(i-1))/ba(i); end % Back substitution step u(R-1,n+1)=ga(R-1); for i=R-2:-1:2 u(i,n+1)=ga(i)-c(i)*u(i+1,n+1)/ba(i); end dt=1.1*dt; end % end of time step loop % % % % Plot Storing plot value of u as up, at every 5 time steps j=time index i=space index for j=5:5:N+1 for i=1:R up(i)=u(i,j); end plot(x,up) hold on end grid title('u[xx]=u[t] 0< x< 0.5; r(u[xx])=u[t] 0.5< x< 1; u(0,t)=1, u(x,0)=0; u[x]=r(u[x]) x=0.5') xlabel('x - ND Space') ylabel('u - ND Temperature') hold off gtext('r=0.01') u(1,t)=1, % Storing times for temp. profiles % These can be saved in a data file or examined in the command file tp=zeros(1,(N-1)/5); i=1; tp(1)=0; for k=5:5:N+1 i=i+1; tp(i)=t(k); end tp tp = Columns 1 through 7 0 0.0029 0.0085 0.0175 0.0320 0.0553 0.0929 0.4079 0.6607 1.0679 1.7238 2.7799 Columns 8 through 14 0.1534 0.2509 Column 15 4.4809 u[x x ]= u[t] 0< x < 0.5; r(u[x x ])= u[t] 0.5< x < 1; u(0,t)= 1, u(1,t)= 1, u(x ,0)= 0; u[x ]= r(u[x ]) x = 0.5 1.4 r= 1 1.2 u - ND Tem perature 1 0.8 0.6 0.4 0.2 0 0 0.1 0.2 0.3 0.4 0.5 0.6 x - ND S pac e 0.7 0.8 0.9 1 u[x x ]= u[t] 0< x < 0.5; r(u[x x ] )= u[t] 0. 5< x < 1; u(0,t)= 1, u(1,t)= 1, u(x ,0)= 0; u[x ]= r(u[x ]) x = 0.5 1. 4 1. 2 r= 0.1 u - ND Tem perature 1 0. 8 0. 6 0. 4 0. 2 0 0 0.1 0.2 0.3 0. 4 0. 5 0.6 0.7 0.8 0. 9 1 x - N D S pac e u[x x ]= u[t] 0< x < 0.5; r(u[x x ] )= u[t] 0. 5< x < 1; u(0,t)= 1, u(1,t)= 1, u(x ,0)= 0; u[x ]= r(u[x ]) x = 0.5 1 0. 9 r= 0.0 1 0. 8 u - ND Tem perature 0. 7 0. 6 0. 5 0. 4 0. 3 0. 2 0. 1 0 0 0.1 0.2 0.3 0. 4 0. 5 0.6 x - N D S pac e 0.7 0.8 0. 9 1 u[x x ]= u[t] 0< x < 0.5; r(u[x x ] )= u[t] 0. 5< x < 1; u(0,t)= 1, u(1,t)= 1, u(x ,0)= 0; u[x ]= r(u[x ]) x = 0.5 1 0. 9 r= 0.001 0. 8 u - ND Tem perature 0. 7 0. 6 0. 5 0. 4 0. 3 0. 2 0. 1 0 0 0.1 0.2 0.3 0. 4 0. 5 0.6 0.7 0.8 0. 9 1 x - N D S pac e u[x x ]= u[t] 0< x < 0.5; r(u[x x ] )= u[t] 0. 5< x < 1; u(0,t)= 1, u(1,t)= 1, u(x ,0)= 0; u[x ]= r(u[x ]) x = 0.5 1. 4 1. 2 u - ND Tem perature 1 r= 0 0. 8 0. 6 0. 4 0. 2 0 0 0.1 0.2 0.3 0. 4 0. 5 0.6 x - N D S pac e 0.7 0.8 0. 9 1 32.14 99.27296 99.38879 99.47967 99.38879 99.27296 25 10 10 10 5 99.15714 99.40126 99.57055 99.40126 99.15714 40 21.87149 13.44564 7.62124 0 dh/dx 1.040287 1.014349 0.931015 dh/dy 0.040287 0.054636 0.068985 dh/dn 1.041067 1.015819 0.933568 θ (radians) 0.038707 0.053811 0.073961 θ (degrees) 2.217773 3.083137 4.237646 16.3372 16.29691 16.22792 20 20 20 40 24.04033 14.28983 7.039322 0 96.07795 100 100 100 96.07795 30 20 12.63401 6.246222 0 85.75874 88.97417 100 88.97417 85.75874 15 10 5.311556 0 1.106512 1.057395 0.778698 1.311258 1.344371 0.62638 1.449779 1.5883 0.066225 0.109272 0.152318 0.138521 0.38245 0.62638 0 0 1.108492 1.063026 0.793455 1.318555 1.397713 0.885835 1.449779 1.5883 0.059779 0.102975 0.193167 0.105249 0.277161 0.785398 0 0 3.425088 5.90002 11.06765 6.030345 15.88014 45 0 0 17.37748 17.31126 17.15894 18.55022 18.4117 17.78532 20 20 69.00866 70.13795 72.56898 70.13795 69.00866 7.5 5 2.5 0 0 0 0 0 0 0 0.052697 0.082316 0.082316 0.052697 0 0 0.072156 0.113101 0.113101 0.072156 0 0 0.052697 0.082316 0.082316 0.052697 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.27719 -0.39124 -0.39124 -0.27719 0 0 -0.35074 -0.50218 -0.50218 -0.35074 0 0 -0.27719 -0.39124 -0.39124 -0.27719 0 0 0 0 0 0 0 1.387741 1.391891 1.409973 c 76.44 52.47 36.06 25.05 19.09 98.55306 99.48858 100 99.48858 98.55306 0.113952 0.168488 0.48078 dcdx -9.588 -8.076 -5.484 -3.394 -2.384 J-diff 19.176 16.152 10.968 6.788 4.768 0.155496 0.793428 6.185917 J-adv 76.44 52.47 36.06 25.05 19.09 0.864874 6.581653 100 J 95.616 68.622 47.028 31.838 23.858 0.951623 6.938975 100 0.958962 6.959813 100 50 50 50 50 50 c0 c1 c2 c3 c4 76.53 63.25 52.28 43.22 35.75 c5 c6 c7 c8 c9 29.61 24.62 20.69 17.88 16.58