[ Pobierz całość w formacie PDF ]
.GT.ABS(XEND)-1.E-8) RETURNCALL RUKU4S(N,X,DX,YTT) ! Solve with full incrementERRM=0.0X1=X+DXDO 30 I=1,NYTT(I)=Y(I)-YTT(I)30 ERRM=MAX(ERRM,ABS(YTT(I)/Y(I)))IF(ERRM.EQ.0) THENDX=5.*DXDXS=DXGO TO 1ELSEERRM=ERRM/ERRORDX=DX/ERRM**0.2DXS=DXIF(ERRM.GT.1.0) THENDO 40 I=1,NYTT(I)=YORI(I)40 Y(I)=YORI(I)GO TO 20ENDIFENDIFDO 50 I=1,N50 Y(I)=Y(I)+YTT(I)/15.! Accounts for truncation errorGO TO 1ENDFigure A.6 A fourth-order Runge-Kutta code with automatic adjustment of step size.The arguments for RUKUST now have the following meanings:N = number of ODE's to be solved, and for which derivatives will be given.© 2000 by CRC Press LLCDXS = an initial value for ∆ x.This value will be decreased or increased, dependingupon what is needed to satisfy the error criterion.In previous Runge-Kutta sub-routines DX was the interval over which the problem was solved.Now DXSnormally will be smaller than this value.Upon returning from this subroutine,DXS is the ∆ x that was found to be satisfactory at the end of the solution, andit can be used as the initial increment in a subsequent call to RUKUST.XBEG = the initial value of the independent variable.XEND = the end value of the independent variable.The difference between XENDand XBEG was called DX in the previous subroutines.ERROR = the error criterion to meet in obtaining the numerical solution.Y = an array of N values that, upon entry to the subroutine, provides the initial con-ditions for the dependent variable.Upon return from the subroutine it is the solu-tion for the dependent variable(s) at x = XEND.YTT = a work array of N values.It is used only to store the solution for the last in-terval, which is then compared with the solution Y in making decisions aboutthe next increment in the independent variable to use in satisfying the errorcriterion.Observe that this subroutine calls RUKU4S three times; the first two times com-plete the solution over the increment ∆ x in two steps of length ∆ x/2 (this solution isstored in array Y), and the third time uses the increment ∆ x, i.e.uses the four sub-stepsin the Runge-Kutta method (and this solution is stored in the work array YTT).Thedifference between these two solutions is used to determine the error ERR (or ERR1);then, based on Eq.A.15, the ∆ x that should supply the desired accuracy is computed.Ifthe accuracy is insufficient, then the solution is repeated, using the computed ∆ x (thestatement GO TO 20 does this).Another test checks whether the solution has proceededto XEND.If not, then the solution proceeds over the next increment with the newlycomputed ∆ x by going back to statement 1.The program is required to end the solutionat XEND, and this is accomplished by adjusting the last ∆ x so it equals the differencebetween the current value of x and XEND.Example Problem A.3The bottom of the tank is defined by the ODE dy/dx = x + y/2, measured in a coordi-nate system that is rotated downward 45o from the horizontal.The bottom is located be-tween two vertical walls that are 2 m apart, and the tank is 5 m long.It contains waterthat is 5 m deep at the left wall.Find the vertical component of force on the bottom, andthe location of its line of action.allW2mallWb = 5m5mWateryy'45°x'dy/dx = x + y/2x© 2000 by CRC Press LLCTo solve this problem, we must first solve the ODE to obtain the shape of the bot-tom and then numerically integrate from this bottom position to the water surface.Theprogram below obtains the solution.We will need to establish the relation between the x-direction and the horizontal direction, denoted by x' in the sketch, to account for the rotat-ed coordinate system.The differential area to integrate can be written dA = hdx' , in whichh is the distance from the bottom of the tank to the water surface; from the sketch thisdistance is h = H + x sin θ - y cos θ, with H = 5 being the water depth at the origin.Also from the sketch x' = x cos θ + y sin θ.Since sin 45o = cos 45o, only sin 45o will be used in the program.While there are alternate approaches, in this solution the ODEwill be solved first to obtain the shape of the bottom, with the results stored in an array,YY; when needed, values from this array will be interpolated.We choose this approachbecause, although the numerical integration is in terms of x' and the corresponding valueof x along the rotated axis will be needed, the only way to determine this x is to usex' /sin 45o = x + y.For any x' , therefore, the table is searched for the entry where x + y is just larger than x' /sin 45o, and then a linear interpolation is used to find x by writingy = yo + (∆ y/∆ x)( x - xo) between the two entries in the table, with subscript o denoting the first entry.Thus x' /sin 45o = yo + (∆ y/∆ x)( x - xo) + x, and the solution is found to be x = ( x' /sin 45o - yo + (∆ y/∆ x) xo)/(∆ y/∆ x + 1).Once x is determined, y is interpolated as y = yo + ( x - xo)/[∆ x( y1 - yo)].The ODE is solved first in the main program, and the function EQUAT performs the interpolations to obtain h so thatSimpson's rule can properly evaluate the area.The solution produces a cross-sectional areaof 9.44 m2, leading to a force per unit length of 92.5 kN/m and a total vertical force FV= 5x92.5 = 462.5 kN.The first moment of the area is determined by changing theEQUAT statement slightly, as the comment statement in the program shows, and theresult of the integration gives xcA = 8.82 m3 so the line of action of this vertical force isat xc = 8.82/9.44 = 0.934 m.TANKODEK [ Pobierz całość w formacie PDF ]