CH-08 (523178), страница 2
Текст из файла (страница 2)
Any component hasa difference exceeding the tdf value will cause flag1 to change from a value of 1 to0 and the Break command in the second Do loop in In[2] to exit and to continuethe iteration. flag1 is created to control the While command which determines whenthe iteration should be terminated. The N[tn,3] instructs the components of {tn} beprinted with 3 significant figures.When tdf is changed to a value of 0.5, Mathematica can again be applied to yieldIn[4]: = t = (Table[0,{n}]; tn = Table[0,{n}]; t[[1]] = 100; tn[[1]] = 100;tdf = 0.5; tm = 0; flag1 = 0;)In[5]: = (While[flag1 = = 0,Do[tn[[j]] = t[[j]] + c1*(t[[j–1]]–2*t[[j]] + t[[j + 1]]),{j,2,nm1}];tn[[n]] = t[[n]] + 2*c1*(t[[nm1]]-t[[n]]);tm = tm + dt;flag1 = 1;Do[If[Abs[tn[[I]]-t[[I]]]>tdf, flag1 = 0; Break,Continue],{i,n}];Do[t[[I]] = tn[[I]],{i,n}]]; Print[“t = “,tm]; Print[N[tn,3]])Out[5]: = t = 49{100., 75.5, 53.2, 34.9, 21.3, 12., 6.24, 3., 1.35, 0.617, 0.413}Notice that Mathematica takes only 49 seconds, one second less than thatrequired for the FORTRAN, QuickBASIC, and MATLAB versions.
The reason isthat Mathematica keeps more significant digits in carrying out all computations.To show more on the effect of changing the tdf value, the following Mathematicaruns are provided:© 2001 by CRC Press LLCIn[6]: = t = (Table[0,{n}]; tn = Table[0,{n}]; t[[1]] = 100; tn[[1]] = 100;tdf = 0.1; tm = 0; flag1 = 0;)In[7]: = (While[flag1 == 0, Do[tn[[j]] = t[[j]] + c1*(t[[j–1]]–2*t[[j]] + t[[j + 1]]),{j,2,nm1}];tn[[n]] = t[[n]] + 2*c1*(t[[nm1]]-t[[n]]);tm = tm + dt; flag1 = 1;Do[If[Abs[tn[[I]]-t[[I]]]>tdf, flag1 = 0; Break,Continue],{i,n}];Do[t[[I]] = tn[[I]],{i,n}]]; Print[“t = “,tm];Print[N[tn,3]])Out[7]: = t = 462{100., 93.9, 88., 82.3, 77.1, 72.5, 68.5, 65.3, 63., 61.6, 61.1}In[8]: = t = (Table[0,{n}]; tn = Table[0,{n}]; t[[1]] = 100; tn[[1]] = 100;tdf = 0.05; tm = 0; flag1 = 0;)In[9]: = (While[flag1 == 0, Do[tn[[j]] = t[[j]] + c1*(t[[j–1]]–2*t[[j]] + t[[j + 1]]),{j,2,nm1}];tn[[n]] = t[[n]] + 2*c1*(t[[nm1]]-t[[n]]);tm = tm + dt; flag1 = 1;Do[If[Abs[tn[[I]]-t[[I]]]>tdf, flag1 = 0; Break,Continue],{i,n}];Do[t[[I]] = tn[[I]],{i,n}]]; Print[“t = “,tm];Print[N[tn,3]])Out[9]: = t = 732{100., 97., 94., 91.2, 88.5, 86.2, 84.2, 82.6, 81.5, 80.8, 80.5}In[10]: = t = (Table[0,{n}]; tn = Table[0,{n}]; t[[1]] = 100; tn[[1]] = 100;tdf = 0.0001; tm = 0; flag1 = 0;In[11]: = (While[flag1 == 0, Do[tn[[j]] = t[[j]] + c1*(t[[j–1]]–2*t[[j]] + t[[j + 1]]),{j,2,nm1}];tn[[n]] = t[[n]] + 2*c1*(t[[nm1]]-t[[n]]);tm = tm + dt; flag1 = 1;Do[If[Abs[tn[[I]]-t[[I]]]>tdf, flag1 = 0; Break,Continue],{i,n}];Do[t[[I]] = tn[[I]],{i,n}]]; Print[“t = “,tm]; Print[N[tn,3]])Out[11]: = t = 3159{100., 100., 100., 100., 100., 100., 100., 100., 100., 100., 100.}For tdf = 0.1 and tdf = 0.05, Mathematica continues to take lesser time than theother version; but when tdf = 0.0001, Mathematica needs two additional seconds.
The© 2001 by CRC Press LLCreason is that seven significant figures are required in the last case, rounding mayhave resulted in earlier termination of the iteration when the FORTRAN, QuickBASIC, and MATLAB versions are employed.8.3 PROGRAM RELAXATN — SOLVING ELLIPTICAL PARTIALDIFFERENTIAL EQUATIONS BY RELAXATION METHODThe program Relaxatn is designed for solving engineering problems which aregoverned by elliptical partial differential equation of the form:∂2 φ ∂2 φ+= F(x, y)∂x 2 ∂y 2(1)where is called the field function and F(x,y) is called forcing function.
When thesteady-state heat conduction of a two-dimensional domain is considered,2 then thefield function becomes the temperature distribution, T(x,y), and the forcing functionbecomes the heat-source function, Q(x,y). If the distribution of T is influenced onlyby the temperatures at the boundary of the domain, then Q(x,y) = 0 and Equation1 which is often called a Poisson equation is reduced to a Laplace equation:∂2 T ∂2 T+=0∂x 2 ∂y 2(2)The second-order, central-difference formulas (read the program DiffTabl) canbe applied to approximate the second derivatives in the above equation at an arbitrarypoint x = xi and y = yj in the domain as:T − 2 Ti, j + Ti +1, j∂2 Tat x i , y j =˙ i −1, j2∂x( ∆x)2andT − 2 Ti, j + Ti, j+1∂2 Tat x i , y j =˙ i, j−12∂y(∆y)2By substituting both of the above equations into Equation 2 and taking equalincrements in both x- and y-directions, the reduced equation is, for x = y,Ti, j =(1T + Ti +1, j + Ti, j−1 + Ti, j+14 i −1, j)(3)The result is expected when the temperature distribution reaches a steady statebecause it states that the temperature at any point should be equal to the average ofits surrounding temperatures.© 2001 by CRC Press LLCFIGURE 3.
A plate, which initially (at time t = 0) had a temperature equal to 0°F throughout,insulated along a portion of its boundary and suddenly heated at its upper left boundary tomaintain a linearly varying temperature Tb.Before we proceed further, it is appropriate at this time to introduce a numericalcase. Suppose that a plate which initially (at time t = 0) has a temperature equal to0°F throughout and is insulated along a portion of its boundary, is suddenly heatedat its upper left boundary to maintain a linearly varying temperature Tb as shown inFigure 3. If this heating process is to be maintained, we are then interested inknowing the temperature distribution, changed from its initial state of uniformlyequal to 0°F if given sufficient time to allow it to reach an equilibrium (steady) state.Numerically, we intend to calculate the temperatures, denoted as a matrix [T], at aselected number of locations.
Therefore, the plate is first divided into a gridwork ofM rows and N columns along the x- and y-directions, respectively, as indicated inFigure 1. The directions of x- and y-axes are so selected for the convenience ofassociating them with the row and column of the temperature matrix [T] which isof order M by N. The values of M and N should be so decided such that theincrements x and y are equal in order to apply Equation 3. To be more specific,let M = 10 and N = 20 and the linear temperature variation along the upper leftboundary Tb be:Ti,1 = 10(i − 1)fori = 1, 2,…, 6(4)Here, Ti,j is to be understood as the temperature at the location (xi,yj).
Equation4 describes the temperature along the left boundary y = y1 but only for xi in therange of i = 1 to i = 6.Since the temperatures at the stations which are on the insulated boundaries ofthe plate are also involved, these unknown temperatures need to be treated differently.By an insulated boundary, it means that there is no heat transfer normal to theboundary. Since the heat flow is depended on the temperature difference across that© 2001 by CRC Press LLCinsulated boundary, mathematically it requires that T/n = 0 there, n being thenormal direction.
At the vertical boundaries x = xM, we have T/n = T/x = 0 sincex is the normal direction. Based on the central difference and considering twoincrements in the x direction, at a yjth station we can have:T− TM −1, j∂Tat x M , y j =˙ M +1, j=0∂x2( ∆x)()orTM +1, j = TM −1, j(5)Since xM + 1 is below the bottom boundary of the plate shown in Figure 3, thereis no need to calculate the temperatures there, however Equation 5 enables thetemperatures at the stations along the bottom boundary of the plate TM,j for j = 1 toN to be averaged.
Returning to Equation 4, we notice that if Equation 5 is substitutedinto it, the resulting equation which relates only to three neighboring temperatures is:TM, j =(12 TM −1, j + TM, j−1 + TM, j+14)for j = 2, 3,…, N − 1(6)Notice that j = 1 and j = N are not covered in Equation 6.
These two casesconcerning the insulated stations at the left and right bottom corners of the platewill be discussed after we address the two vertical, insulated boundaries y = y1 andy = yN.For the boundaries y = y1, ∂T/∂n becomes ∂T/∂y. Again, we can apply the centraldifference for double y increments to obtain:T − Ti,2∂Tat (x1, y1 ) =˙ i,0=0∂y2( ∆y)orTi,0 = Ti,2(7)Thus, the modified Equation 3 for the left insulated boundary is:Ti,1 =(1T + Ti +1,1 + 2 Ti,24 i −1,1)fori = 7, 8, 9In a similar manner, we can derive for the right insulated boundary y = yNT− Ti,N −1∂Tat (x i , y N ) =˙ i,N +1=0∂y2( ∆y)or© 2001 by CRC Press LLC(8)Ti,N +1 = Ti,N −1(9)andTi,N =(1T+ Ti +1,N + 2 Ti,N −14 i −1,N)for i = 8, 9(10)Having derived Equations 6, 8, and 10, it is easy to deduce the two specialequation for the corner insulated stations to be:TM,1 =(1T + TM −1,12 M ,2)(11)andTM,N =(1T+ TM −1,N2 M,N −1)(12)We have derived all equations needed for averaging the temperature at any stationof interest including those at the insulated boundaries by utilizing those at itsneighboring stations.
It suggests that a continuous upgrading process can be developed which assumes that the neighboring temperatures are known. This so-calledrelaxation method starts with an initial assumed distribution of temperature [T(0)]and continues to use Equations 3 and 6 to 12 until the differences at all locationsare small enough. Mathematically, the process terminates when:MNi =1j=1∑∑ T( k +1)i, j− Ti(,kj ) < ε(13)where is a prescribed tolerance of accuracy and k = 0,1,2,… is the number ofsweeps in upgrading the temperature distribution.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
