CH-06 (523176), страница 5
Текст из файла (страница 5)
This solution is only for the selected value of N. Nshould be continuously increased to test if the maximum Z value would be affected.The program OdeBvpFD has been developed in both QuickBASIC and FORTRAN versions based on the above described procedure for solving the boundaryproblems governed by ordinary differential equations.
It is designed for the generalcase where the coefficient matrix [C] and the right-hand side vector {R} of thematrix equation derived from the finite-difference approximation, [C]{Z} = {R}, are© 2001 by CRC Press LLCFIGURE 12. As a second example, we consider the simply supported beam of length 3L whichhas a solid section at its middle portion and hollow sections at both end portions as shownto be defined by the user. In addition, the user needs to specify the number of station,N, the stepsize, R, and the location of the left boundary, Ri. Here, Z and R arereferred as the dependent and independent variables, respectively using the membrane problem only for the convenience of explanation; user could have othernotations as dependent and independent variables.As a second example, we consider the simply supported beam of length 3Lwhich has a solid section at its middle portion and hollow sections at both endportions as shown in Figure 12.
It is of interest to know how the beam will deflectunder its own weight. This is the case of a beam subjected to uniformly distributedloads. For simplicity, let us assume that the distributed loads to be wm and we, inN/cm2, for the middle and end portions, respectively. To determine the deflectiony(x) for 0<x<3L, we need the following relevant equations which are available in astandard textbook on mechanics of materials7:dy= θ,dxdθ d 2 y M== ,dx dx 2 EIdM= V,dxdV=wdx(12–15)where , M, and V are the slope, bending moment, and shearing force, respectively.E is the Young’s modulus and I is the moment of inertia of the beam. The distributedload w is considered as positive when it is applied upward in the direction of thepositive y-axis and the moment is considered as positive if it causes the beam tobend concave up.
The shearing forces are considered as positive if they are relatedto M according to Equation 14.By successive substitutions, Equations 12–15 can yield an equation which relatesthe deflection y directly to the distributed load w as:d2 d2y EI =wdx 2 dx 2 © 2001 by CRC Press LLC(16)This is a fourth-order ordinary differential equation and needs fourth supplementary initial, boundary, or mixed conditions in order to completely solve for thedeflection y(x).
For the simply supported beam subjected to its own weight, the fourboundary conditions are:y=0andM = EId2y=0dx 2atx=0andx=L(17)If EI is not a function of x, the central-difference approximation for Equation16 at x = xk is:y − 4 y k −1 + 6y k − 4 y k +1 + y k +2 wd 4y=at x k k −2at x k4dxEI( ∆x)4(18)If N in-between stations are selected for determination of the deflections there,then the two boundaries are at k = 0 and k = N + 1. In view of the subscripts k–2and k + 2 in Equation 18, k can thus only take the values 2 through N–1. For solvingy1 through yN, we hence need two more equations which are the two momentconditions in Equation 17.
At x0 = 0, the second-order, forward-difference formulacan be applied to giveM = EIy − 2 y1 + y 0d2yat x 0 2EI at x 0 = 02dx( ∆x)2Since y = 0 at the left support x = x0 = 0, the above equation yields−2 y1 + y 2 = 0(19)Similarly, at the right support x = xN + 1 = 3L, the second-order, backwarddifference formula can be applied to giveM = EIy − 2 y n + y N −1d2yat x N +1 N +1EI at x N +1 = 02dx( ∆x)2Since y = 0 at the right support x = xN + 1 = 3L, the above equation yieldsy N −1 − 2 y N = 0(20)Meanwhile, the boundary condition y = 0 at the left support x = x0 = 0 can besubstituted into Equation 18 for k = 2 to obtain:−4 y1 + 6y 2 − 4 y3 + y 4 = ( ∆x)© 2001 by CRC Press LLC4wat x 2EI(21)Similarly, the boundary condition y = 0 at the right support x = xN + 1 = 3L canbe substituted into Equation 18 for k = N–1 to obtain:y N −3 − 4 y N −2 + 6y N −1 − 4 y N = ( ∆x)4wat x N −1EI(22)Equation 18 now should have the modified form of, for k = 3,4,…,N–2:y k −2 − 4 y k −1 + 6y k − 4 y k +1 + y k +2 = ( ∆x)4wat x kEI(23)In matrix form, Equations 19–23 can be written as [C]{Y} = {R} where {Y} =[y1 y2 • • • yN]T and the elements of the coefficient matrix [C] denoted as ci,j are tobe calculated using the formulas:c1,1 = c N,N = −2,ci,i = 6,ci,i −1 = ci,i +1 = −4,ci,i −2 = ci,i +2 = 1,ci, j = 0c1,2 = c N,N −1 = 1for i = 2, 3,…, N − 1for i = 2, 3,…, N − 1for i = 3, 4,…, N − 2elsewhere(24)For the right-side vector {R} in the matrix equation [C]{Y} = {R}, its elementsdenoted as ri are to be calculated using the formula:r1 = rN = 0andri = ( ∆x)4wat x i for i = 2, 3,…, N − 1EI(25)where w(xi) means the distributed load at x = xi.
Consider the beam shown inFigure 12 with dead loads wm and we. That is:w(x i ) = w e′for0 ≤ xi < Lw(x i ) = w m′forL < xi < 2Lw(x i ) = (w m + w e ) 2forxi = Land2 L < x i ≤ 3L(26)andxi = 2LNotice that the average value of wm and we at the junctions of hollow and solidportions of the beam. Expressions for the EI product can be given similar to thosefor w. In a sample application of the QuickBASIC version of the program OdeBvpFD, we will give two subprogram functions which are prepared for the beamproblem shown in Figure 12 based on Equations 25 and 26.An alternative approach for solving the simply supported beam by finite-difference approximation of Equation 13 and using the two boundary conditions y = 0 at© 2001 by CRC Press LLCx = 0 and x = 3L is given as a homework problem.
This approach requires themoment equation M(x) be derived prior to the solution of the matrix equation[C]{Y} = {R}.FORTRAN VERSIONSample ApplicationConsider the membrane problem shown in Figure 1. Let the inner radius Ri beequal to 3 in, the stepsize or radial increment R be equal to 9.5 in, and the numberof stations N between the two boundaries be selected as equal to 11 (that is, theother boundary is at R = Ri + (N + 1) R = 3 + 12x0.5 = Ro = 9 in.) If the tensionT is equal to 100 lbs/in and the pressure p is equal to 5 lbs/in2, a function subprogramRI can then be accordingly prepared as listed in the program OdeBvpFD.
Aninteractive run of this program has resulted in a display on screen as follows:© 2001 by CRC Press LLCQUICKBASIC VERSION© 2001 by CRC Press LLCSample ApplicationsConsider the same membrane problem as for the sample application using theFORTRAN version of the program OdeBvpFD. The QuickBASIC version has aCOMMON SHARED statement allowing N to be shared in the subprograms. Also,this version has been expanded to solving up to 119 stations.
It gives an interactiverun as follows:The simply supported beam under its own weight is considered as a secondexample. For wm = we = –2 and a uniform EI value equal to 1 (that is when the beamhas a uniform cross section without the hollow end portions), L = 100, x = 30,and N = 9 stations between the supports, the subprograms CIJ and RI are preparedaccording to Equations 24 to 26 as follows:© 2001 by CRC Press LLCIt should be noted that the arguments X and DX for FUNCTION CIJ are notactually involved in any of the statements therein.
They are kept so that the statementsin the program OdeBvpFD involving CIJ can remain as general as possible toaccommodate other problems which may need these arguments as linkages. Theinteractive application of the above two FUNCTION subprograms is demonstratedbelow.The analytical solution of this beam problem can be easily obtained.
From anytextbook on mechanics of materials (see footnote *), the maximum deflection for abeam of length 3L and subjected to a uniformly distributed load w is ymax =5w(3L)4/384EI. Since w = –2, L = 100, and EI = 1 are used in the above illustrativerun, ymax is equal to 5x(–2)x(300)4/384 = 2.109x108. The computed result of y at x7(for N = 9, the 5th station is the mid-length of the beam) by the program OdeBvpFDis equal to –1.21500x108 which has an error about 42%. A homework problem isgiven for readers to exercise different N and x values, which will show that for Nequal to 19, 29, 59, 99, and 119, the respective ymax values are –1.63285x108,–1.78535x108, — 1.93665x108, –1.92011x108, and –1.83x108.
Notice that the accuracy continue to improve until the roundoff errors begin to affect the solution whenN becomes large, for such cases the double precision should be used in solving thematrix equation [C]{Y} = {R} by the Gaussian elimination method.MATLAB APPLICATIONSThe two sample problems discussed in FORTRAN and QuickBASIC versionof the program OdeBvpFD can be executed interactively by MATLAB with thecommands entered from keyboard as follows:Membrane Problem© 2001 by CRC Press LLC© 2001 by CRC Press LLCBeam Deflection Problem© 2001 by CRC Press LLCMATHEMATICA APPLICATIONSFor solving the membrane problem, Mathematica can be applied as follows:In[1]: = Ns = 11; DX = 0.5; c = Table[0,{Ns},{Ns}];R = Table[–5*DX^2/100,{Ns}];Print[“R = “,R]Out[1] = R = {–0.0125, –0.0125, –0.0125, –0.0125, –0.0125, –0.0125,–0.0125, –0.0125, –0.0125, –0.0125, –0.0125}In[2]: = (Do[Do[x = 3 + i*DX; If[i = = j, c[[i,j]] = –2;,If[i = = j–1, c[[i,j]] = 1DX/2/x;,If[i = = j + 1, c[[i,j]] = 1 + DX/2/x;, Continue]]],{i,Ns}], {j,Ns}]); Print[“Matrix c = “,c]Out[2] = Matrix c = {{–2, 0.928571, 0, 0, 0, 0, 0, 0, 0, 0, 0},{1.0625,–2, 0.9375, 0, 0, 0, 0, 0, 0, 0, 0},{0, 1.05556, –2, 0.944444, 0, 0, 0, 0, 0, 0, 0},{0, 0, 1.05, –2, 0.95, 0, 0, 0, 0, 0, 0},{0, 0, 0, 1.04545, –2, 0.954545, 0, 0, 0, 0, 0},{0, 0, 0, 0, 1.04167, –2, 0.958333, 0, 0, 0, 0},{0, 0, 0, 0, 0, 1.03846, –2, 0.961538, 0, 0, 0},{0, 0, 0, 0, 0, 0, 1.03571, –2, 0.964286, 0, 0},{0, 0, 0, 0, 0, 0, 0, 1.03333, –2, 0.966667, 0},{0, 0, 0, 0, 0, 0, 0, 0, 1.03125, –2, 0.96875},{0, 0, 0, 0, 0, 0, 0, 0, 0, 1.02941, –2}}In[3]: = V = Inverse[c].ROut[3] = {0.0533741, 0.101498, 0.142705, 0.175525, 0.198641, 0.210864,0.211107, 0.198368, 0.171723, 0.13031, 0.0733212}These results are in agreement with those obtained by the MATLAB application.The loaded beam problem also can be treated in a similar manner as follows:In[1]: = c = Table[0,{9},{9}]; c[[1,1]] = –2; c[[1,2]] = 1; c[[9,9]] = –2;c[[9,8]] = 1;In[2]: = (Do[Do[If[i = = j,c[[i,j]] = 6;, If[(j = = i–1)||(j = = i + 1), c[[i,j]] = –4;,If[(j = = i + 2)||(j = = i–2), c[[i,j]] = 1;, Continue]]],{i,2,8}],{j,9}]); Print[“Matrix c = “,c]Out[2] = Matrix c = {{–2, 1, 0, 0, 0, 0, 0, 0, 0}, {–4, 6, –4, 1, 0, 0, 0, 0, 0},{1, –4, 6, –4, 1, 0, 0, 0, 0}, {0, 1, –4, 6, –4, 1, 0, 0, 0},{0, 0, 1, –4, 6, –4, 1, 0, 0}, {0, 0, 0, 1, –4, 6, –4, 1, 0},© 2001 by CRC Press LLC{0, 0, 0, 0, 1, –4, 6, –4, 1}, {0, 0, 0, 0, 0, 1, –4, 6, –4},{0, 0, 0, 0, 0, 0, 0, 1, –2}}In[3]: = R = Table[0,{9}]; WE = –2; WM = –2; EIE = 1; EIM = 1; L = 100;DX = 30;In[4]: = (Do[X = i*DX; If[X<L, R[[i]] = DX^4*WE/EIE;,If[X = = L, R[[i]] = DX^4*(WE + WM)/(EIE + EIM);,If[X>L, If[X<2*L, R[[i]] = DX^4*WM/EIM;,If[X = = 2*L, R[[i]] = DX^4*(WE + WM)/(EIE + EIM);,IF[X>2*L, R[[i]] = DX^4*WE/EIE;,Continue]]]]]],{i,2,8}]); Print[“R = “,R]Out[4] = R = {0, –1620000, –1620000, –1620000, –1620000, –1620000,–1620000, –1620000, 0}In[5]: = V = Inverse[c].ROut[5] = {–34020000, –68040000, –96390000, –115020000, –121500000,–115020000, –96390000, –68040000, –34020000}Again, the results are in agreement with those obtained by the MATLABapplication.6.5 PROBLEMSRUNGEKUT1.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
