Chapter 31 The Equation to Be Solved is

CHAPTER 31 31.1 The equation to be solved is d 2T dx 2 = −20 Assume a solution of the form T = ax2 + bx + c which can be differentiated twice to give T” = 2a. Substituting this result into the differential equation gives a = −10. The boundary conditions can be used to evaluate the remaining coefficients. For the first condition at x = 0, 50 = −10(0) 2 + b(0) + c or c = 50. Similarly, for the second condition. 100 = −10(10) 2 + b(10) + 50 which can be solved for b = 105. Therefore, the final solution is T = −10 x 2 + 105x + 50 The results are plotted in Fig. 31.5. 300 200 100 0 0 5 10 31.2 The heat source term in the first row of Eq. (31.26) can be evaluated by substituting Eq. (31.3) and integrating to give 2 .5 ∫0 20 2.5 − x dx = 25 2.5 Similarly, Eq. (31.4) can be substituted into the heat source term of the second row of Eq. (31.26), which can also be integrated to yield 2 .5 ∫0 20 x−0 dx = 25 2.5 These results along with the other parameter values can be substituted into Eq. (31.26) to give 0.4T1 − 0.4T2 = − and dT ( x ) + 25 dx 1 −0.4T1 + 0.4T2 = dT ( x ) + 25 dx 2 31.3 In a manner similar to Fig. 31.7, the equations can be assembled for the total system, −0.4 0.8  0.4  −0.4      −0.4 −0.4 0.8 −0.4 −0.4 0.8 −0.4   T1   T  2      T3  =   −0.4 T4   0.4    T5   − dT ( x 1 ) / dx + 25   50     50     50     dT ( x1 ) / dx + 25   The unknown end temperatures can be substituted to give 1       −0.4 0.8 −0.4 −0.4 0.8   dT ( x1 ) / dx   5     70  T2          T3  =  50     90  T4     −1  − dT ( x 5 ) / dx     −15  −0.4 −0.4 0.8 −0.4 These equations can be solved for  dT ( x1 ) / dx    T2     T3  =   T4    − dT ( x5 ) / dx   105  250      325  275     −95  The solution, along with the analytical solution (dashed line) is shown below: 300 200 100 0 0 5 31.4 0=D R= D ∫ x2 x1 d 2c dx 2 −U dc − kc dx ~ d 2~ c dc ~ −U − kc 2 dx dx  d 2~ c dc~ ~  N dx −U − kc D  i 2 dx dx   10 D x2 d 2~ c x1 dx ∫ −U x2 ∫x 1 −k ∫ x2 N i ( x ) dx (1) ~ dc N ( x )dx dx i (2) 2 ~ N ( x ) dx c i x1 (3) Term (1): D ∫ c − c2   dc − ( x1 ) − 1  d c  dx x 2 − x1   N i ( x ) dx = D  2 c 2 − c1  dc dx   ( x2 ) −  x 2 − x1   dx  2~ x2 x1 Term (2): ∫ dc~ N ( x ) dx = dx i x2 x1 x2 ∫x N i ( x ) dx = 1 ∴ 1 ∫ x2 x1 c2 − c1 N ( x) dx x 2 − x1 i x 2 − x1 2 dc~ c − c1 N ( x ) dx = 2 dx i 2 x2 ∫x −U ∫ x2 x1 c2 ~  dc N i ( x )dx = −U  c dx  2  − c1   2 − c1   2  Term (3): −k x2 ∫x 1 k ( x 2 − x1 ) c~N i ( x ) dx = − 2  c1     c2  Total element equation [(1) + (2) + (3)] a 11 a  11 a 11  c1   b1   =  a 11 c 2  b2  where a11 = D U k − + ( x 2 − x1 ) x 2 − x1 2 2 a 22 = − D U k + + ( x 2 − x1 ) x2 − x1 2 2 a12 = −D U + x 2 − x1 2 a 21 = −D U − x 2 − x1 2 b1 = − D dc ( x1 ) dx b2 = D dc ( x2 ) dx 31.5 First we can develop an analytical function for comparison. Substituting parameters gives 15 . × 10 8 2 d u dx 2 = 50 Assume a solution of the form u = ax 2 + bx + c This can be differentiated twice to yield d2u/dx2 = 2a. This can be substituted into the ODE, which can then be solved for a = 1.6667×10−7. The boundary conditions can then be used to evaluate the remaining coefficients. At the left side, u(0) = 0 and 0 = 16667 . × 10 −7 ( 0) 2 + b(0) + c and therefore, c = 0. At the right side of the bar, u(30) = 0 and 0 = 16667 . × 10 −7 ( 30) 2 + b( 30) and therefore, b = −5×10−6, and the solution is u = 1.6667 × 10 −7 x 2 − 5 × 10 −6 x which can be displayed as 0 6 12 18 24 30 0E+00 -2E-05 -4E-05 The element equation can be written as  du   Ac E  1 −1  u1  − dx ( x1 )   +     = Ac E  du x 2 − x1  −1 1  u2   (x 2 )    dx   x2 ∫x x ∫x 1 1 2 P( x ) N 1 ( x ) dx   P ( x ) N 2 ( x ) dx   The distributed load can be evaluated as 6 ∫0−50 6−x dx = −150 6 6 ∫0−50 Thus, the final element equation is x−0 dx = −150 6  2.5 × 10 7  7  −2.5 × 10  du  −2.5 × 10 7   u1  − dx ( x 1 )  −150     = Ac E  du +  −150  2.5 × 10 7  u 2   ( x2 )    dx  Assembly yields 1.5 × 108         −2.5 × 10 7 7 5 × 10 −2.5 × 10 7 −2.5 × 10 5 × 10 7 7 −2.5 × 10 7 −2.5 × 10 7 7 5 × 10 −2.5 × 10 7 −2.5 × 10 5 × 10 7 7 −2.5 × 10 7 which can be solved for  du    ( x1 )  −150  dx       u1  −300    u1  −300   =    u1  −300    u  −300  1     7  du −1.5 × 10   ( x 2 ) −150   dx   du  −6  dx ( x1 )   −5 × 10     u  −2.4 × 10 −5  1     u1   −3.6 × 10 −5   = −5   u1   −3.6 × 10   u1   −2.4 × 10 −5   du   −6    ( x 2 )   5 × 10  dx  These results, along with the analytical solution (dashed line) are displayed below: 0 6 12 18 24 30 0E+00 -2E-05 -4E-05 31.6 Option Explicit Sub FErod() Dim Dim Dim Dim Dim Dim Dim Dim ns As Integer, ii As Integer, i As Integer, j As Integer k As Integer, m As Integer x(5) As Single, st(2, 2) As Single, c As Single s(2, 2) As Single, a(5, 5) As Single, b(5) As Single, d(5) As Single Te(5) As Single, ff As Single e(5) As Single, f(5) As Single, g(5) As Single, r(5) As Single dum1 As Single, dum2 As Single dTeLeft As Single, dTeRight As Single 'set parameters ns = 4 x(1) = 0 x(2) = 2.5 x(3) = 5 x(4) = 7.5 x(5) = 10 Te(1) = 40 Te(5) = 200 ff = 10 'construct system matrix st(1, 1) = 1: st(1, 2) = -1: st(2, 1) = -1: st(2, 2) = 1 For ii = 1 To ns c = 1 / (x(ii + 1) - x(ii)) For i = 1 To 2 For j = 1 To 2 s(i, j) = c * st(i, j) Next j Next i For i = 1 To 2 k = ii - 1 + i For j = 1 To 2 m = ii - 1 + j a(k, m) = a(k, m) + s(i, j) Next j b(k) = b(k) + ff * ((x(ii + 1) - x(ii)) - (x(ii + 1) - x(ii)) / 2) Next i Next ii 'determine impact of uniform source and boundary conditions Call Mmult(a(), Te(), d(), ns + 1, ns + 1, 1) For i = 1 To ns + 1 b(i) = b(i) - d(i) Next i a(1, 1) = 1 a(2, 1) = 0 a(ns + 1, ns + 1) = -1 a(ns, ns + 1) = 0 'Transform square matrix into tridiagonal form f(1) = a(1, 1) g(1) = a(1, 2) r(1) = b(1) For i = 2 To ns e(i) = a(i, i - 1) f(i) = a(i, i) g(i) = a(i, i + 1) r(i) = b(i) Next i e(ns + 1) = a(ns + 1, ns) f(ns + 1) = a(ns + 1, ns + 1) r(ns + 1) = b(ns + 1) 'Tridiagonal solver dum1 = Te(1) dum2 = Te(ns + 1) Call Tridiag(e, f, g, r, ns + 1, Te()) dTeLeft = Te(1) dTeRight = Te(ns + 1) Te(1) = dum1 Te(ns + 1) = dum2 'output results Range("a3").Select ActiveCell.Value = "dTe(x = " & x(0) & ")/dx = " ActiveCell.Offset(0, 1).Select ActiveCell.Value = dTeLeft ActiveCell.Offset(1, -1).Select ActiveCell.Value = "dTe(x = " & x(ns + 1) & ")/dx = " ActiveCell.Offset(0, 1).Select ActiveCell.Value = dTeRight ActiveCell.Offset(3, -1).Select For i = 1 To ns + 1 ActiveCell.Value = x(i) ActiveCell.Offset(0, 1).Select ActiveCell.Value = Te(i) ActiveCell.Offset(1, -1).Select Next i Range("b3").Select End Sub Sub Mmult(a, b, c, m, n, l) Dim i As Integer, j As Integer, k As Integer Dim sum As Single For i = 1 To n sum = 0 For k = 1 To m sum = sum + a(i, k) * b(k) Next k c(i) = sum Next i End Sub Sub Tridiag(e, f, g, r, n, x) Dim k As Integer 'decompose For k = 2 To n e(k) = e(k) / f(k - 1) f(k) = f(k) - e(k) * g(k - 1) Next k 'substitute For k = 2 To n r(k) = r(k) - e(k) * r(k - 1) Next k x(n) = r(n) / f(n) For k = n - 1 To 1 Step -1 x(k) = (r(k) - g(k) * x(k + 1)) / f(k) Next k End Sub The output is 31.7 After setting up the original spreadsheet, the following modifications would be made to insulate the right edge and add the sink: Cell I1: Set to 100 Cell I2: =(I1+2*H2+I3)/4; This formula would then be copied to cells I3:I8. Cell I9: =(I8+H9)/2 Cell C7: =(C6+D7+C8+B7-110)/4 The resulting spreadsheet is displayed below: A 87.5 75 75 75 75 75 75 75 75 1 2 3 4 5 6 7 8 9 B 100 89.6 86.4 85.0 82.2 75.6 66.8 70.3 72.0 C 100 96.9 96.2 96.3 93.2 78.4 46.1 67.5 72.2 D 100 101.7 105.2 110.8 115.7 98.8 81.8 81.4 82.0 E F 100 104.9 112.1 126.0 160.1 119.4 100.8 94.3 92.8 100 105.7 112.4 120.9 129.4 117.9 107.7 102.3 100.7 G 100 105.4 110.8 115.9 118.7 114.9 109.9 106.5 105.3 H 100 105.1 109.5 113.1 114.6 113.2 110.5 108.4 107.6 I 102.8 105.6 109.2 112.2 113.5 112.6 110.6 108.9 108.2 Corresponding contour plots can be generated as S9 S8 S7 S6 S5 S4 S3 S2 1 2 3 4 5 6 7 8 9 S1 180 160 140 120 100 80 60 40 20 0 S9 S7 S5 S3 S1 1 2 3 4 5 6 7 8 9 31.8 The results of the preceding problem (31.8) can be saved as a tab-delimited text file (in our case, we called the file prob3108.txt). The following commands can then be used to load this file into MATLAB, as well as to generate the contour plot along with heat flow vectors. >> >> >> >> load prob3108.txt [px,py]=gradient(prob3108); cs=contour(prob3108);clabel(cs);hold on quiver(-px,-py);hold off 9 80 8 7 6 120 140 160 5 4 100 60 3 2 1 1 2 3 4 5 6 7 8 9 31.9 The scheme for implementing this calculation on Excel is shown below: 1 2 3 4 5 6 A 87.5 75 75 75 75 62.5 B 100 C 100 D 100 E 100 F 100 G 100 H 100 I 100 J 100 50 50 50 50 50 50 50 50 50 K 62.5 25 25 25 25 37.5 The simple Laplace equation applies directly to the blank white cells (recall Fig. 31.14). However, for the shaded cells impacted by the heat sink, a heat balance equation must be written. For example, for cell E3, the control volume approach can be developed as 0 = −k' E 3 − D3 F 3 − D3 E3 − E4 E 2 − E3 ∆y∆z + k ' ∆y∆z − k ' ∆ x∆z + k ' ∆x∆ z − 100 ∆x∆ y ∆x ∆x ∆y ∆y Collecting and canceling terms yields 0 = −4 E 3 + D3 + F 3 + E 4 + D 3 − 100 ∆x∆y ∆zk ' Substituting the length dimensions and the coefficient of thermal conductivity gives, E3 = D 3 + F 3 + E 4 + D3 − 160 4 The result is depicted below, along with the corresponding contour plots. 1 2 3 4 5 6 A 87.5 75 75 75 75 62.5 B 100 74.0 59.0 53.1 53.3 50 C 100 62.0 33.7 25.3 34.9 50 D E F G 100 100 100 100 40.2 10.1 -1.5 7.9 -11.2 -98.2 -123.8 -101.8 -20.7 -108.0 -133.7 -111.6 11.2 -19.5 -31.3 -21.8 50 50 50 50 H 100 34.7 -19.8 -29.2 5.7 50 I 100 50.8 17.3 8.8 23.8 50 J 100 51.3 29.2 23.4 30.5 50 K 62.5 25 25 25 25 37.5 S6 S5 2 3 4 5 6 7 8 9 10 S1 11 S3 S1 11 9 7 1 S5 5 S2 3 S3 100 50 0 -50 -100 -150 1 S4 31.10 The results of the preceding problem (31.10) can be saved as a tab-delimited text file (in our case, we called the file prob3110.txt). The following commands can then be used to load this file into MATLAB, as well as to generate the contour plot along with heat flow vectors. >> >> >> >> load prob3110.txt [px,py]=gradient(prob3110); cs=contour(prob3110);clabel(cs);hold on quiver(-px,-py);hold off 6 5.5 5 4.5 -50 -100 4 3.5 3 2.5 2 0 1.5 1 31.11 2 50 4 6 8 10 Program Plate Use IMSL Implicit None Integer::ncval, nx, nxtabl, ny, nytabl Parameter (ncval=11, nx=33, nxtabl=5, ny=33, nytabl=5) Integer::i, ibcty(4), iorder, j, nout Real::ax,ay,brhs,bx,by,coefu,prhs,u(nx,ny),utabl,x,xdata(nx),y,ydata(ny) External brhs, prhs ax = 0 bx = 40 ay = 0 by = 40 ibcty(1) = 1 ibcty(2) = 2 ibcty(3) = 1 ibcty(4) = 1 coefu = 0 iorder = 4 Call FPS2H(prhs, brhs, coefu, nx, ny, ax, bx, ay, by, ibcty, iorder, u, nx) Do i=1, nx xdata(i) = ax + (bx - ax) * Float(i - 1) / Float(nx - 1) End Do Do j=1, ny ydata(j) = ay + (by - ay) * Float(j - 1) / Float(ny - 1) End Do Call UMACH(2, nout) Write (nout,'(8X,A,11X,A,11X,A)') 'X', 'Y', 'U' Do j=1, nytabl Do i=1, nxtabl x = ax + (bx - ax) * Float(i - 1) / Float(nxtabl - 1) y = ay + (by - ay) * Float(j - 1) / Float(nytabl - 1) utabl = QD2VL(x,y,nx,xdata,ny,ydata,u,nx,.FALSE.) Write (nout,'(4F12.4)') x, y, utabl End Do End Do End Program Function prhs(x, y) Implicit None Real::prhs, x, y prhs = 0 End Function Real Function brhs(iside, x, y) Implicit None Integer::iside Real::x , y If (iside == 1) Then brhs = 50 ElseIf (iside == 2) Then brhs = 0 ElseIf (iside == 3) Then brhs = 75 Else brhs = 100 End If End Function Output: 0.0000 10.0000 20.0000 30.0000 40.0000 0.0000 10.0000 20.0000 30.0000 40.0000 0.0000 10.0000 20.0000 30.0000 40.0000 0.0000 10.0000 20.0000 30.0000 40.0000 0.0000 10.0000 20.0000 30.0000 40.0000 Press any key to 0.0000 0.0000 0.0000 0.0000 0.0000 10.0000 10.0000 10.0000 10.0000 10.0000 20.0000 20.0000 20.0000 20.0000 20.0000 30.0000 30.0000 30.0000 30.0000 30.0000 40.0000 40.0000 40.0000 40.0000 40.0000 continue 75.0000 71.6339 66.6152 59.1933 50.0000 75.0000 72.5423 67.9412 60.1914 50.0000 75.0000 75.8115 72.6947 64.0001 50.0000 75.0000 83.5385 83.0789 74.3008 50.0000 87.5000 100.0000 100.0000 100.0000 75.0000 31.12 Other node equations are derived similarly 31.13 Other node equations are derived similarly