Chapter 29 After 15 Iterations of the Liebmann Method, The Result is

CHAPTER 29 29.1 First iteration: 7.500000 9.750000 55.425000 Error: 100.000000 100.000000 100.000000 Second iteration: 9.600000 26.137500 68.068500 Error: 21.875000 62.697270 18.574660 • • • 2.250000 3.600000 62.707500 15.675000 20.782500 85.047000 100.000000 100.000000 100.000000 100.000000 100.000000 100.000000 8.212501 34.632000 88.782750 20.563500 52.916250 85.500310 72.602740 89.604990 29.369720 23.772710 60.725670 5.301830E-01 Seventh iteration: 25.013610 28.806340 46.216590 56.257030 78.575310 93.082440 Error: 2.954020E-01 2.531316E-02 2.267362E-02 2.082395E-02 2.165254E-03 3.590016E-03 33.932440 56.921290 87.501180 1.679560E-02 1.041445E-02 1.743838E-03 29.2 The fluxes for Prob. 29.1 can be calculated as qx= -9.325527E-02 -7.657973E-01 -1.668020 qy= -1.132306 -1.312262 -2.542694 qn= 1.136140 1.519367 3.040984 theta= -94.708180 -120.266600 -123.265100 -2.185114E-01 -2.622653E-01 -2.186839E-01 -5.192447E-01 1.532972E-01 1.055520 -1.378297 -1.574765 -2.296703 -1.394572 -1.312434 -2.280428 1.395511 1.596454 2.307091 1.488101 1.321357 2.512862 -99.008540 -99.455400 -95.439100 -110.421900 -83.337820 -65.162450 29.3 The plate is redrawn below 100 oC 75 oC 50 oC 0 oC After 15 iterations of the Liebmann method, the result is 0 50 50 50 50 50 50 50 0 100 73.6954 62.3814 55.9987 51.1222 46.0804 39.2206 27.1773 0 100 82.3973 69.8296 60.4898 52.4078 43.9764 33.6217 19.4897 0 100 86.06219 74.0507 63.72554 54.04625 43.79945 31.80514 17.16646 0 100 87.7991 76.58772 66.32058 56.25934 45.37425 32.62971 17.3681 0 100 88.54443 78.18341 68.71677 59.3027 48.80563 35.95756 19.66293 0 100 88.19118 78.8869 71.06672 63.42793 54.57569 42.71618 25.31308 0 100 85.32617 78.10995 73.23512 68.75568 63.33804 54.995 38.86852 0 0 75 75 75 75 75 75 75 0 0 0.0070 0.0079 0.0007 0.0067 0.0077 0.0630 0.2194 0 0 0.0114 0.0120 0.0097 0.0164 0.0400 0.1192 0.2925 0 0 0.0120 0.0109 0.0241 0.0542 0.1005 0.2343 0.7119 0 0 0 0 0 0 0 0 0 0 with percent approximate errors of 0 0.0030 0.0050 0.0062 0.0076 0.0106 0.0149 0.0136 0 0 0 0 0 0 0 0 0 0 0 0.0040 0.0062 0.0067 0.0066 0.0079 0.0099 0.0013 0 0 0.0043 0.0057 0.0036 0.0020 0.0033 0.0119 0.0302 0 0 0.0049 0.0055 0.0007 0.0106 0.0074 0.0327 0.1259 0 29.4 The solution is identical to Prob. 29.3, except that now the top edge must be modeled. This means that the nodes along the top edge are simulated with equations of the form 4Ti , j − Ti −1, j − Ti +1, j − 2 Ti , j −1 = 0 The resulting simulation (after 14 iterations) yields 50 50 50 50 50 50 50 50 0 50.38683 50.17211 49.51849 48.31607 46.33449 43.09381 37.46764 26.36368 0 51.16385 50.76425 49.56564 47.39348 43.91569 38.56608 30.4051 17.98153 0 52.6796 52.15054 50.58556 47.78093 43.37764 36.8614 27.61994 15.18654 0 with percent approximate errors of 55.17802 54.58934 52.86931 49.79691 44.99165 37.93565 28.08718 15.20479 0 58.7692 58.20129 56.56024 53.61405 48.94264 41.91332 31.71478 17.63115 0 63.41846 62.96008 61.64839 59.2695 55.38806 49.21507 39.39338 23.73251 0 68.9398 68.67918 67.93951 66.58047 64.29121 60.37012 53.1291 38.00928 0 75 75 75 75 75 75 75 75 0 0 0 0 0 0 0 0 0 0 0.0584 0.0722 0.0854 0.0933 0.0930 0.0873 0.0913 0.1131 0 0.1318 0.1603 0.1883 0.2121 0.2300 0.2469 0.2827 0.3612 0 0.2034 0.2473 0.2937 0.3441 0.3913 0.4299 0.4995 0.7054 0 0.2606 0.3173 0.3788 0.4464 0.5097 0.5474 0.5852 0.9164 0 0.2828 0.3424 0.4077 0.4754 0.5328 0.5611 0.5525 0.7958 0 0.2493 0.2983 0.3438 0.3972 0.4468 0.4237 0.3157 0.5085 0 0.1529 0.1862 0.2096 0.2247 0.2605 0.2747 0.0477 0.6345 0 0 0 0 0 0 0 0 0 0 29.5 The solution is identical to Examples 29.1 and 29.3, except that now heat balances must be developed for the three interior nodes on the bottom edge. For example, using the control-volume approach, node 1,0 can be modeled as −0.49 (5) T10 − T00 10 + 0.49 (5) T20 − T10 10 + 0.49(10) T11 − T10 10 − 2(10) = 0 4T10 − T00 − T20 − 2 T11 = −81.63265 The resulting simulation yields (with a stopping criterion of 1% and a relaxation coefficient of 1.5) 87.5 75 75 75 75 100 79.91669 66.88654 52.26597 27.12079 100 77.76103 60.34068 40.84576 10.54741 100 70.67812 55.39378 40.26148 14.83802 The fluxes for the computed nodes can be computed as qx -0.06765 0.359153 0.836779 1.579088 0.226345 0.281573 0.29411 0.300928 0.680145 0.253347 -0.22428 -0.96659 qy -0.81128 -0.67744 -0.97426 -1.23211 -0.97165 -0.90442 -1.21994 -1.48462 -1.09285 -0.74521 -0.99362 -1.24575 qn 0.814095 0.766759 1.284283 2.002904 0.997668 0.947241 1.254887 1.514811 1.287216 0.787095 1.018614 1.576764 θ (radians) -1.65398 -1.08331 -0.86117 -1.34193 -1.26898 -1.33422 -1.0141 -1.24309 -1.7928 75 50 50 50 50 -0.66259 -1.37081 -2.23067 θ (degrees) -94.7663 -62.0692 -49.3412 -37.9638 -76.8869 -72.7072 -76.4454 -78.5416 -58.1036 -71.2236 -102.72 -127.808 29.6 The solution is identical to Example 29.5 and 29.3, except that now heat balances must be developed for the interior nodes at the lower left and the upper right edges. The balances for nodes 1,1 and 3,3 can be written as − 4 T11 + 0.8453T21 + 0.8453T12 = −1154701 . ( T01 + T10 ) − 4 T33 + 0.8453T32 + 0.8453T23 = −1154701 . (T34 + T43 ) Using the appropriate boundary conditions, simple Laplacians can be used for the remaining interior nodes. The resulting simulation yields 75 100 100 100 50 75 86.02317 94.09269 100 50 63.97683 75 86.02317 100 50 55.90731 63.97683 75 100 50 50 50 75 29.7 The nodes to be simulated are 0,3 1,3 2,3 3,3 0,2 1,2 2,2 3,2 0,1 1,1 2,1 3,1 0,0 1,0 2,0 3,0 Simple Laplacians are used for all interior nodes. Balances for the edges must take insulation into account. Fo example, node 1,0 is modeled as 4T1,0 − T0,0 − T2 ,0 − 2T1,1 = 0 The corner node, 0,0 would be modeled as 4T0, 0 − 2 T1,0 − 2T0,1 = 0 The resulting set of equations can be solved for 0 11.94853 15.625 16.36029 16.36029 12.5 25 37.5 16.08456 22.79412 30.14706 17.09559 19.94485 22.79412 16.72794 17.09559 16.08456 16.36029 15.625 11.94853 50 37.5 25 12.5 0 The fluxes can be computed as Jx -0.6125 -0.20267 -0.07206 -0.01801 -5.6E-13 -0.6125 -0.26572 -0.10584 -0.01801 0.018015 -0.6125 -0.34453 -0.13961 0.015763 0.108088 -0.6125 -0.36029 -0.12385 0.112592 0.382812 -0.6125 -0.36029 -0.10809 0.175643 0.585478 Jy 0.585478 0.382812 0.108088 0.018015 0 0.175643 0.112592 0.015763 -0.01801 -0.01801 -0.10809 -0.12385 -0.13961 -0.10584 -0.07206 -0.36029 -0.36029 -0.34453 -0.26572 -0.20267 -0.6125 -0.6125 -0.6125 -0.6125 -0.6125 Jn 0.847314 0.43315 0.129906 0.025477 5.63E-13 0.637187 0.288587 0.107004 0.025477 0.025477 0.621964 0.366116 0.197444 0.107004 0.129906 0.710611 0.509533 0.366116 0.288587 0.43315 0.866206 0.710611 0.621964 0.637187 0.847314 θ (radians) 2.378747 2.057696 2.158799 2.356194 3.141593 2.862322 2.740799 2.993743 -2.35619 -0.7854 -2.96692 -2.7965 -2.35619 -1.42295 -0.588 -2.60987 -2.35619 -1.91589 -1.17 -0.4869 -2.35619 -2.10252 -1.74547 -1.29153 -0.80795 θ (degrees) 136.2922 117.8973 123.6901 135 180 163.999 157.0362 171.5289 -135 -45 -169.992 -160.228 -135 -81.5289 -33.6901 -149.534 -135 -109.772 -67.0362 -27.8973 -135 -120.466 -100.008 -73.999 -46.2922 29.8 Node 0,3: There are two approaches for modeling this node. One would be to consider it a Dirichlet node and not model it at all (i.e., set it’s temperature at 50oC). The second alternative is to use a heat balance to model it as shown here (0,3) (1,3) (0,2) 0 = 0.5(15)(1) T1, 3 − T0,3 40 − 0.5( 20)(1) T0,3 − T0,2 30 + 0.01(20 )(1)(10 − T0,3 ) −0.29752 T1,3 + 4T0 ,3 − 0.52893T0,2 = 3.17355 Node 2,3: 0 = −0.5(15)(1) 4T2 ,3 T2 ,3 − T1,3 + 0.5(15)(1) T3,3 − T2 ,3 40 20 − 0.70588 T1,3 − 1.41176 T3,3 − 188235 . T2 ,2 − 0.5( 30)(1) T2, 3 − T2, 2 30 + 0.01( 30)(1)(10 − T2 ,3 ) Node 2,2: (2,3) (2,2) (3,2) (1,2) (2,1) 0 = −0.5( 22.5)(1) T2 ,2 − T1,2 40 + 0.5(22.5)(1) T3,2 − T2 ,2 20 − 0.5( 30)(1) T2 ,2 − T2,1 15 T2 ,3 − T2 ,2 + 0.5( 30)(1) + 10π ( 7.5) 2 30 4T2 ,2 − 0.48T1,2 − 0.96 T3,2 − 170667 . T2,1 − 0.85333T2 ,3 = 3015.93 Node 5,3: (5,3) (4,3) (5,2) 0 = −0.5(15)(1) T5,3 − T4 ,3 20 − 0.5(10)(1) T5,3 − T5,2 30 + 0.01(10)(1)(10 − T5,3 ) 4T5,3 − 2.33766 T4,3 − 1.03896T5,2 = 6.23377 29.9 Node 0,0: (0 ,1 ) (1 ,0 ) (0 ,0 ) 0 = 0.01( 7.5)( 2)( 20 − T0, 0 ) − 0.7( 7.5)( 2 ) T1,0 − T0,0 40 + 0.7( 20)( 2) T0,1 − T0,0 15 4T0, 0 − 0.46069 T1, 0 − 3.27605T0,1 = 5.26508 Node 1,1: (1,2) (1,1) (2,1) (0,1) (1,0) 0 = −0.7 ( 22.5)( 2) T1,1 − T0,1 + 0.5( 22.5)(2 ) T2 ,1 − T1,1 − 0.7 (20 )(2 ) T1,1 − T1, 0 40 20 15 T1,1 − T1, 0 T1,2 − T1,1 T1,2 − T1,1 − 0.5(10 )(2 ) + 0.7( 20)( 2) + 0.5(10)( 2) 15 30 30 4T1,1 − 0.78755T2,1 − 1.77389 T1, 0 − 0.88694 T1,2 − 0.55142 T0,1 = 0 Node 2,1: (2,2) (2,1) (1,1) (3,1) (2,0) 0 = −0.5( 22.5)( 2 ) T2 ,1 − T1,1 20 + 0.5( 22.5)( 2) T3,1 − T2 ,1 20 − 0.5( 20)(2 ) − 0.5( 20)( 2) T2 ,1 − T2 ,0 T2, 2 − T2,1 30 15 + 10( 22.5)(20) 4T2 ,1 − 1.05882 T1,1 − 1.05882 T3,1 − 1.2549 T2,0 − 0.62745T2 ,2 = 4235.29 29.10 The control volume is drawn as in 0, j+1 1, j 0, j-1 A flux balance around the node can be written as (note ∆x = ∆y = h) − kh∆z T1, j − T0, j h + k (h / 2 ) ∆z T0, j − T0, j −1 h Collecting and cancelling terms gices T0, j − T0, j −1 − T0, j +1 − 2T1, j = 0 − k (h / 2 ) ∆z T1, j − T0, j h =0 29.11 A setup similar to Fig. 29.11, but with θ > 45o can be drawn as in 2 θ 3 8 4 1 ∆y 7 6 5 ∆x The normal derivative at node 3 can be approximated by the gradient between nodes 1 and 7, T −T ∂T = 1 7 ∂η 3 L17 When θ is greater than 45o as shown, the distance from node 5 to 7 is ∆y cotθ, and linear interpolation can be used to estimate T7 = T5 + (T6 − T5 ) ∆y cot θ ∆x The length L17 is equal to ∆y/sinθ. This length, along with the approximation for T7 can be substituted into the gradient equation to give  ∆y  ∂ T T1 =    sin θ  ∂ η − T6 3 ∆ y cot θ ∆y cot θ   − T5  1 −   ∆x ∆x  29.12 The following Fortran-90 program implements Liebmann’s method with relaxation. PROGRAM liebmann IMPLICIT NONE INTEGER :: nx,ny,l,i,j REAL :: T(0:5,0:5),ea(0:5,0:5),Told(0:5,0:5) REAL :: qy(0:5,0:5),qx(0:5,0:5),qn(0:5,0:5),th(0:5,0:5) REAL :: Trit,Tlef,Ttop,Tbot,lam,emax,es,pi REAL :: k,x,y,dx,dy nx=4 ny=4 pi=4.*atan(1.) x=40. y=40. k=0.49 lam=1.2 es=1. dx=x/nx dy=y/ny Tbot=0. Tlef=25. Trit=50. Ttop=150. DO i=1,nx-1 T(i,0)=Tbot END DO DO i=1,nx-1 T(i,ny)=Ttop END DO DO j=1,ny-1 T(0,j)=Tlef END DO DO j=1,ny-1 T(nx,j)=Trit END DO l=0 DO l=l+1 emax=0. DO j = 1,ny-1 DO i = 1,nx-1 Told(i,j)=T(i,j) T(i,j)=(T(i+1,j)+T(i-1,j)+T(i,j+1)+T(i,j-1))/4 T(i,j)=lam*T(i,j)+(1-lam)*Told(i,j) ea(i,j)=abs((T(i,j)-Told(i,j))/T(i,j))*100. if(ea(i,j).GT.emax) emax=ea(i,j) END DO END DO PRINT *, 'iteration = ',l DO j = 1,ny-1 PRINT *, (T(i,j),i=1,nx-1) END DO PRINT * DO j = 1,ny-1 PRINT *, (ea(i,j),i=1,nx-1) END DO IF (emax.LE.es) EXIT END DO DO j = 1,ny-1 DO i = 1,nx-1 qy(i,j)=-k*(T(i,j+1)-T(i,j-1))/2/dy qx(i,j)=-k*(T(i+1,j)-T(i-1,j))/2/dx qn(i,j)=sqrt(qy(i,j)**2+qx(i,j)**2) th(i,j)=atan2(qy(i,j),qx(i,j))*180./pi END DO END DO PRINT *,'qx=' DO j = 1,ny-1 PRINT *, (qx(i,j),i=1,nx-1) END DO PRINT *,'qy=' DO j = 1,ny-1 PRINT *, (qy(i,j),i=1,nx-1) END DO PRINT *,'qn=' DO j = 1,ny-1 PRINT *, (qn(i,j),i=1,nx-1) END DO PRINT *,'theta=' DO j = 1,ny-1 PRINT *, (th(i,j),i=1,nx-1) END DO END When the program is run, the result of the last iteration is: iteration = 42.81303 63.17175 78.57594 0.5462000 1.1274090E-02 3.1769749E-02 qx= 1.022510 0.4589829 -2.7459882E-02 qy= -1.547708 -0.8761914 -0.9022922 qn= 1.854974 0.9891292 0.9027100 theta= -56.54881 -62.35271 -91.74317 Press any key to 6 33.26489 56.26600 76.12081 33.93646 52.46138 69.64268 0.1074174 2.0983342E-02 3.6572997E-02 2.4864437E-02 4.8064217E-02 2.4659829E-02 0.2174759 0.2624041 0.2188648 -0.4100102 0.1535171 0.6399599 -1.378517 -1.049970 -1.071483 -1.285304 -0.8748025 -1.164696 1.395566 1.082263 1.093608 1.349116 0.8881705 1.328934 -81.03486 -75.96829 -78.45538 continue -107.6926 -80.04664 -61.21275 29.13 When the program is run, the result of the last iteration is: iteration = 25.01361 46.21659 78.57531 7 28.80634 56.25703 93.08244 33.93244 56.92129 87.50118 0.2954020 2.5313158E-02 1.6795604E-02 2.2673620E-02 2.0823948E-02 1.0414450E-02 2.1652540E-03 3.5900162E-03 1.7438381E-03 qx= -9.3255267E-02 -0.2185114 -0.5192447 -0.7657973 -0.2622653 0.1532972 -1.668020 -0.2186839 1.055520 qy= -1.132306 -1.378297 -1.394572 -1.312262 -1.574765 -1.312434 -2.542694 -2.296703 -2.280428 qn= 1.136140 1.395511 1.488101 1.519367 1.596454 1.321357 3.040984 2.307091 2.512862 theta= -94.70818 -99.00854 -110.4219 -120.2666 -99.45540 -83.33782 -123.2651 -95.43910 -65.16245 Press any key to continue 29.14 When the program is run, the result of the last iteration is: iteration = 19 38.490 24.764 19.044 54.430 41.832 34.955 62.710 53.570 47.660 68.165 62.478 58.219 72.761 70.301 67.841 16.783 31.682 44.291 55.234 65.489 16.696 31.041 42.924 53.215 63.051 19.176 33.110 43.388 51.848 60.034 27.028 38.976 45.799 50.854 55.780 77.795 78.373 77.594 76.027 73.595 85.175 87.944 88.261 87.530 85.843 69.522 62.233 82.249 73.624 This data can be imported into Excel and the following contour plot created: 40 35 30 25 20 15 80-100 60-80 40-60 20-40 0-20 10 5 0 5 10 15 20 25 30 35 0 40 29.15 29.16