CH-08 (523178), страница 4
Текст из файла (страница 4)
When the total temperature differences, d, is limited to eps = 100 degrees, the t values obtained after twosweeps are slightly different from those obtained by the other versions, this is againbecause Mathematica keeps more significant digits in all computation steps thanthose in the FORTRAN, QuickBASIC, and MATLAB programs. By changing theeps value from 100 degrees to 1 degree, Mathematica also takes 137 sweeps toconverge as in the FORTRAN, QuickBASIC, and MATLAB versions:In[5]: = t = Table[0,{10},{20}]; eps = 1; nr = 0; d = eps + 1;In[6]: = Do[t[[i,1]] = (I–1)*10,{i,2,6}];In[7]: = Print[“Sweep #”,nr]; Round[N[t,2]]Out[7] = Sweep #137{{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0},{10, 8, 7, 5, 4, 3, 3, 2, 2, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0},{20, 16, 13, 10, 8, 7, 5, 4, 3, 3, 2, 2, 2, 1, 1, 1, 1, 0, 0, 0},{30, 24, 19, 15, 12, 9, 8, 6, 5, 4, 3, 3, 2, 2, 1, 1, 1, 1, 0, 0},{40, 30, 23, 18, 15, 12, 10, 8, 6, 5, 4, 4, 3, 2, 2, 1, 1, 1, 0, 0},{50, 34, 26, 20, 16, 13, 11, 9, 7, 6, 5, 4, 3, 3, 2, 2, 1, 1, 0, 0},{35, 30, 25, 21, 17, 15, 12, 10, 8, 7, 6, 5, 4, 3, 3, 2, 2, 1, 1, 0},{29, 27, 24, 21, 18, 15, 13, 11, 9, 7, 6, 5, 4, 3, 3, 2, 2, 1, 1, 1},{26, 26, 23, 21, 18, 16, 13, 11, 9, 8, 6, 5, 4, 4, 3, 2, 2, 1, 1, 1},{26, 25, 23, 21, 18, 16, 13, 11, 9, 8, 7, 5, 4, 4, 3, 2, 2, 2, 1, 1}}The above results all agree with those obtained earlier.WARPING ANALYSISOF ATWISTED BARWITHNONCIRCULAR CROSS SECTIONAs another example of applying the relaxation method for engineering analysis,consider the case of a long bar of uniform rectangular cross section twisted by thetwo equal torques (T) at its ends.
The cross section of the twisted bar becomeswarped as shown in Figure 6. If z-axis is assigned to the longitudinal direction ofthe bar, to find the amount of warping at any point (x,y) of the cross sectionedsurface, denoted as W(x,y), the relaxation method can again be employed because© 2001 by CRC Press LLCFIGURE 6.
In the case of a bar of uniform rectangular cross section twisted by the twoequal torques (T) at its ends, the cross section of the twisted bar becomes warped as shown.W(x,y) is governed by the Laplace equation.3 Due to anti-symmetry of W(x,y), thatis W(–x,y) = W(x,–y) = –W(x,y), only one-fourth of the cross section needs to beanalyzed. Let us consider a rectangular region 0≤x≤a and 0≤y≤b.
It is obvious thatthe anti-symmetry leads to the conditions W(x,0) = 0 and W(0,y) = 0 along two ofthe four linear boundaries. To derive the boundary conditions along the right sidex = a and the upper side y = b, we have to utilize the relationships between thewarping function W(x,y) and shear stresses xz and yz which are: ∂Wτ xz = Gθ− y ∂xand ∂Wτ yz = Gθ+ x ∂y(19,20)where G and are the shear rigidity and twisting angle of the bar, respectively.Along the boundaries x = a and y = b, the shear stress should be tangential to thelateral surface of the twisted bar requires: dx ∂W dy ∂W− y−=0+ x ∂x ds ∂y ds(21)where s is the variable changed along the boundary. Along the upper boundary y =b, dy/ds = 0 and dx/ds = 1 and along the right boundary x = a, dx/ds = 0 and dy/ds =1. Consequently, Equation 19 reduces to:∂W= −x∂yfor y = band∂W=yxfor x = a(22,23)To apply the relaxation method for solving W(x,y) which satisfies the Laplaceequation for 0<x<a and 0<y<b and the boundary conditions W(0,y) = W(x,0) = 0© 2001 by CRC Press LLCand Equations 22 and 23, let us partition the rectangular region 0≤x≤a and 0≤y≤binto a network of MxN using the increments x = a/(M–1) and y = b/(N–1).
Infinite-difference forms, Equations 22 and 23 yield, respectively:W(x i , y = y N ≡ b) = W(x i , b − ∆y) − x i ∆yfor i = 2, 3,…, M − 1(24)j = 2, 3,…, N − 1(25)and()()W x = x M ≡ a, y j = W a − ∆x, y j + y j∆xforwhere xi = (i–1)x and yj = (j–1)y. Equations 24 and 25 can be combined to yielda finite-difference formula for relaxing the corner point (xM,yN). That is:[W(x M ≡ a, y N ≡ b) = .5 W(a, y N −1 ) + W(x M −1, b) + b∆x − a∆y](26)Program Relaxatn.m has been modified to develop a new m file Warping.mfor the warping analysis which uses input values for a, b, M, N, and which is neededin Equation 13 for termination of the relaxation. MATLAB statements for sampleapplication of Warping.m are listed below along with the m-file itself.
The meshplot of the warping surface W(x,y) obtained after 131 sweeps of relaxation byinitially assuming W(x,y) = 0 throughout the entire region and an value equal to100 is shown in Figure 7.FIGURE 7. The result for = 100 when a = 30 and b = 20 means that each of the 600gridpoints is allowed to have, on average, a difference equal to 1/6 during two consecutiverelaxations.© 2001 by CRC Press LLCNotice that the variable Nrelax keeps track how many sweeps have been performed during the relaxation procedure, it is an output argument like W which canbe printed if desired.
The exit flag, FlagExit, is initially set equal to zero and changedto a value of unity when the total error, SumOfDs, is less than and allowing therelaxation to be terminated.Figure 7 is the result for = 100 when a = 30 and b = 20 which means that eachof the 600 gridpoints is allowed to have, on average, a difference equal to 1/6 duringtwo consecutive relaxations. The W values are found to be in the range of –154 and24.8 for = 100.
When is set equal to 1 which means that each gridpoint is allowedto have, on the average, a difference of 1/600, the mesh plot of W is shown inFigure 8. The W values are now all equal to or less than zero.© 2001 by CRC Press LLCFIGURE 8.FIGURE 9.It is noteworthy that if the cross section of the twisted rod is square, that is whena = b, there is no warping along the x = y line as illustrated by the mesh plot shown inFigure 9. This case is left as a homework problem for the readers to work out the details.© 2001 by CRC Press LLCFIGURE 10. The problem of a tightened string.8.4 PROGRAM WAVEPDE — NUMERICAL SOLUTION OF WAVEPROBLEMS GOVERNED BY HYPERBOLIC PARTIALDIFFERENTIAL EQUATIONSProgram WavePDE is designed for numerical solution of the wave problems governed by the hyperbolic partial differential equation.1 Consider the problem of atightened string shown in Figure 10. The lateral displacement u satisfies the equation:∂2 u∂2 u= a2 22∂t∂x(1)where x is the variable along the longitudinal direction of the string and a is thewave velocity related to the tension T in the string and the mass m of the string bythe equation:a2 =Tm(2)Together with the governing equation of u, Equation 1, there are also so-calledinitial conditions prescribed which may be expressed as:u(t = 0, x) = f (x)and∂u(t = 0, x)= v0 (x)∂t(3,4)f(x) in Equation 3 describe the initial position while v0(x) describes the initialvelocity of the tightened string.
And, also there are so-called boundary conditionswhich in the case of a string at both ends are:u(t, x = x1 = 0) = 0© 2001 by CRC Press LLCandu( t , x = x N = L ) = 0(5,6)If the string is made of a single material, T/m would then be equal to a constant.Analytical solution can be found for this simple case. For the general case that thestring may be composed of a number of different materials and the mass m is thena function of the spatial variable x. The more complicated the variation of theseproperties in x and t, the more likely no analytical solution possible and the problemcan only be solved numerically.The finite-difference approximation of Equation 1 can be achieved by applyingthe central difference for the second derivatives for both with respect to the spacevariable x and the time variable t.
If we consider the displacement u only a finitenumber of stations, say N, defined with a spatial increment x and using a timeincrement of t, then specially, for t at ti and x at xj, we can have:∂2 u u i −1, j − 2 u i, j − u i +1, j=∂t 2( ∆t )2and∂2 u u i, j−1 − 2 u i, j + u i, j+1=∂x 2( ∆x)2Substituting the above expressions into Equation 1, we obtain:2a∆t u i +1, j = − u i =1, j + 2 u i, j + u − 2 u i, j + u i, j+1 ∆x i, j−1()(7)Equation 7 is to be applied for j = 2 through j = N–1.
Initially for t = t1 = 0,Equation 7 can be applied for i = 2, then ui–1,j≡ u(t1,xj) term is simply f(xj) and canbe calculated using Equation 3, and the ui,j–1, ui,j, and ui,j + 1 terms can be calculatedusing the forward-difference approximations of Equation 4 which are, for k =2,3,…,N–1u(t 2 , x k ) = u(t1, x k ) + v0 (x k )∆t = f (x k ) + v0 (x k )∆t(8)Since both f(x) and v0(x) are prescribed, use of Equation 7 enables all u to becalculated at t = t3 and at all in-between stations xj, for j = 2,3,…,N–1.
When uvalues at t = t2 and t = t3 are completely known, Equation 7 can again be utilized tocompute u at t = t4 for i = 4 and so forth. Because the string is tightened, we expectthe magnitudes of u values to continuously decrease in time. That is, for a specifiedtolerance, we may demand that the computation be terminated whenmax. u i, j < εj=2~N-1© 2001 by CRC Press LLC(9)It should be noted that the obtained u values are only for the selected incrementst and x.
Whether or not the results will change if either increment is reduced,need be tested.Equation 7 relates the displacement at xj at t = ti + 1 to the displacements at thesame location xj at two previous instants ti and ti–1, and also those of its left andright neighboring points at one time increment earlier, t = ti. This is an approximationwhich ignores the influences of the displacements of its left and right neighboringpoints at the present time t = ti + 1, namely ui + 1,j–1 and ui + 1,j + 1. These influences canbe taken into consideration if the central-difference approximation for the curvatureterm in the x domain is applied for t = ti + 1 instead of at t = ti in derivation of Equation7. Reader is urged to work Problem 5 to find the effect of this change.A computer program called WavePDE has been developed for generating thedeflection of the string at N stations for any specified time increment t until thedeflection are almost all equal to zero throughout.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
