- Email: [email protected]
Solution bounds for varying geometry beams
Journal of Sound and Vibration (1976) 44(3), 43 l-448
SOLUTION
BOUNDS
FOR
T. J. MCDANIEL
VARYING AND
GEOMETRY
BEAMS
V. R. MURTHY
Department of Aerospace Engineering and Engineering Research Institute, Iowa State University, Ames, Iowa 50010, U.S.A. (Received 7 April 1975, and in revised form 15 July 1975)
The theory of differential and integral inequalities is applied to obtain upper and lower bounds to the transfer matrix for beams with varying geometry. Various techniques of generating and refining these bounds are investigated. Numerical results indicate that these bounds can be refined to produce numerical agreement of the upper and the lower bound to a given number of significant digits. Proceeding from bounds on the transfer matrix elements a theory is developed for determining upper and lower bounds on the natural frequencies and mode shapes and on the solution state vector for static loading of such beams. This procedure is then extended to the analysis of multispan beams with varying geometry. Numerical results are presented for various configurations.
1. INTRODUCTION
A recent survey article by McDaniel and Henderson [l] indicates the variety of multispan structures that have been analyzed by the transfer matrix method. The theoretical and numerical aspects of transfer matrix analysis have been developed by a number of researchers, including Pestel and Leckie [2], Lin [3, 41, Mercer [5], Mead [6, 71, Henderson [8, 93, and McDaniel [lo, 111, for application to the dynamics of stiffened beams, plates, and shells. Structures with varying geometry, i.e., varying thickness, varying curvature, etc., can be considered in the transfer matrix analysis by approximating the varying geometry structure by constant geometry segments or by using numerical techniques such as Runge-Kutta, predictor-corrector, or the integrating matrix procedure [12] to obtain numerical transfer matrices. Due to the difficulty of solving variable coefficient differential equations associated with the varying geometry structure, few exact solutions are available for comparison with these techniques. An exact solution is known for simply supported cylindrical shells with varying curvature only for the case of exponentially varying curvature [ 131. The use of small integration steps and double precision arithmetic is standard practice for generating numerical transfer matrices. Even so, it is difficult to determine the accuracy of the numerical method and to detect cases where the accuracy of the numerical method is insufficient. It is important that the numerical procedures have sufficient accuracy but not excessive accuracy, as this increases significantly the cost of computation. If limited accuracy computations are used to determine the individual transfer matrices for a single or multispan structure, how does this affect the final solution? In the following sections this problem will be explored by developing upper and lower bounds for the individual transfer matrices. These bounds will then be used to determine bounds on the final solution and thus to determine the effect of limited accuracy intermediate computations on the final results. 431
T. J. MCDANIEL
432
AND V. R. MURTHY
2. DIFFERENTIAL
EQUALITIES
The theory of differential inequalities can be used to establish bounds on the solution of differential and integral equations and thus can provide answers to the above questions. The texts of Walter [14], Protter and Weinberger [I 51, and Lakshmikantham and Leela [ 161 summarize the main results and provide numerous references. For engineering application, the capability of establishing both upper and lower bounds through the theory of inequalities is quite important. It allows one to bracket unknown exact solutions by rough upper and lower bounds if this is sufficient and by refined upper and lower bounds which agree to a prescribed number of digits when required. Unlike energy methods which insure that complete continuity of the displacement field will produce lower bounds on the static influence coefficients [17] or that complete continuity of the stress field in a structure provides upper bounds to the static influence coefficients, the method of differential and integral inequalities provides both upper and lower bounds to the transfer matrix elements for both the static or harmonically vibrating structure. This suggests that it may be possible to construct upper and lower bounds on natural frequencies and on the corresponding mode shapes whereas the Rayleigh-Ritz and the method of Bazley, Fox, and Stadter [18] yield bounds on only natural frequencies. The primary engineering application of the theory of differential and integral inequalities is to construct bounds on linear and non-linear initial value problems. Since the transfer matrix approach to structural analysis converts the two point linear boundary value problem into an initial value problem, the inequality theorems apply directly. The inequality theorems have not been developed sufficiently to allow the construction of bounds for non-linear, two point boundary value problems except for special classes of problems. Even though extensive development of the theoretical aspects of differential inequalities can be found, few engineering application papers are accompanied by such bounds. Some applications, however, can be cited such as those of Nickel [19] and Adams [20], who have explored boundary layer flow equations using this approach. John [21] has developed error estimates for a non-linear plate analysis, and Cannon [22] has found bounds to a heat conduction problem. In general, high precision in the computed results was not demonstrated. However, in a recent article Bell and Appl [23] developed a computational procedure based on differential inequalities which provided quite accurate bounds. Upper and lower bounds for the response of a non-linear single degree of freedom system agreed to several digits over the range of computation for cases where the system was excited by short duration excitation or provided with given initial conditions. The differential inequalities used in the present beam analysis are of the following form. Consider a linear initial value problem :
;
u-(x)1 = [4-~)lmx)l,
O
(1)
with initial conditions
P-m = [II, and with the elements
of [A(x)], i.e., aij(x), satisfying L7ij
If it is assumed
that aij is piecewise continuous
>
0.
and bounded,
(3) then existence
and uniqueness
BOUNDS FOR VARYING GEOMETRY BEAMS
433
of the solution is assured. Then if [F(X)] and [i(x)] satisfy the following differential inequalities,
wmx)l,
m’(o)12 VI,
(4)
G Lw1m1,
[WV G VI,
(5)
& VP)1 2
ml
$ it follows that
ml
2
.,
Lwl
2
mx)l.
(6)
A proof of this result can be found in Walter’s book [14]. In fact, the theorems presented in this text apply to the case where the elements of [A] may be negative and to non-linear differential equations of the initial value type. Since the elements of [A] are non-negative in the present static or dynamic beam analysis, it can be shown that the solution [T(x)] of equations (1) and (2) has non-negative elements. This can be shown by representing the solution [T(x)] in terms of the matrizant of [A(x)]. This fact is quite useful in the construction of the inequalities (4) and (5). 3. BOUNDS FOR THE STATIC DEFLECTION OF BEAMS One can apply the differential inequalities discussed in the last section to determine bounds on the behavior of beams with variable geometry. Since an exact solution for the transfer matrix for a beam with exponential width variation and constant thickness shown in Figure 1
Figure 1. Beam with exponentialvariation of width. b(x) = b,,exp(ax/l); t(x) = rO; E(x) = Eo. easily obtained, this solution can be compared, as a first example, with bounds on the transfer matrix elements. The exact solution for the transfer matrix is given by
is
‘1
x
t,,
t,, -
1
l23
?24
0
001
x
000
1
(7)
where -ax/l}/a - l/a) I/EZc f&) = (x + t14(x) = (1 + exp {-cl.x/1})Z2x/EZo ~1~+ (exp {-ax/f} - 1) 2f3/EZoa3, [23(x)
=
(1
t24(x)
=
(l/a
-
exp -
{-a//})
(l/a
El0 = Eobo t,3,‘12.
+
&f.?z,, x)
exp
tl,
{-Rx//})
&??&
a,
434
T. J. MCDANIEL
AND V. R. MURTHY
For the cantilever beam with a tip load P shown in Figure analysis that the state vector at station x is given by
WCC, II/(Xl
1, one finds via the transfer
.t,31-t,,
III 1 M,,,
=
matrix
t2J - t24 l-x
V(X)
P.
(8)
-1
where W, $I,Mand Vare usual notation for the deflection, slope, moment, and shear contained in the state vector for a beam. Equation (8) is obtained by substituting boundary conditions W(0) =$(O) = M(I) = 0 and V(I) = -P in the transfer matrix equation (7), solving for the state vector at x = 0, and transmitting that state vector to station x. The initial state vector can be obtained in this particular case from statics and is precisely known. Now, let us apply differential inequalities to obtain bounds on the state vector as a function of x. A comparison of bounds can be made with the exact solution only in special examples such as the present one. In most cases one will be able to compare only the upper and lower bounds. To develop bounds on the transfer matrix, one proceeds from the governing differential equation
Recall that the transfer
with the initial which satisfies
condition
matrix also satisfies the same diherential
equation,
[T(O)] = [I]. Let us choose an approximate
& [Rx)] 2 [A(x)][Qx>]
in
i.e.,
transfer
matrix
[T(x)]
o
and [p(O)] = [r(O)] = [I]. A ccording to the previous section it follows that [F(X)] is an upper bound to [T(x)] for the specified range of x. It is easily shown for the case u ) 0 that the solution of equation (10) with CL= 0 yields a [p(x)] which satisfies the above conditions and is given by t^13= x2/2EI,,
(11) iz4 = x2/2EI,. This upper bound transfer matrix is clearly generated by choosing the smallest El(x) in the given interval. It follows that a lower bound transfer matrix [T(x)] is obtained by choosing the largest value of El(x), i.e., EI,exp{crx,/l}, in the given interval and solving the resulting
BOUNDS
FOR VARYING
GEOMETRY
435
BEAMS
constant coefficient differential equation. This lower bound satisfies the required conditions
and [f(O)] = [T(O)] = [I]. The advantage of the present method of choosing the upper and lower bound transfer matrices is that the required inequalities are easily shown to be satisfied. The disadvantage is that if El(x) varies rapidly along the structure, the range over which [F] and [f] are acceptable bounds is considerably reduced: i.e., x0 is small. A coarse bound on {Z(x)} for the present example is obtained by letting x,, = 1 so that [fWl and PI” x )I are bounds over the entire length of the beam. The resulting boundary state vectors (p(x)} and {k(x)} are given by 1x2/2EZ0- x3j6EZo for
u > 0.
(12)
{i(x)} is obtained in this case by replacing El0 by EZoexp{ctx,/l} in equation (12). One can easily verify that {z} and {i} are upper and lower bounds to {Z} in equation (8). If c1is negative then the upper and lower bounds of equation (12) switch roles. The bounds obtained are coarse since the transfer matrix bounds from which they are derived are not particularly close bounds. Improved upper and lower bounds can be constructed by subdividing the beam into N segments and finding upper and lower bounds to the transfer matrix in each interval. Since for the static or dynamic analysis of beams the elements of a transfer mtarix are non-negative, i.e., tij 2 0, it follows that the upper and lower transfer matrices overj segments are given by PI = rf+,lv+~-Il~~~ rf+21m,
[?I = [f,][i;_,].
. . [7’23[TJ.
(124
As in the previous example, the maximum and minimum values of q in each interval were used to construct the upper and lower bound transfer matrices for that interval. The beam of Figure 1 was divided into ten segments, and the resulting bounds on deflection and slope are plotted as solid lines in Figures 2 and 3. Dividing the beam into fifty segments produced the bounds on deflection and slope given by dashed lines in the same figures. In the latter case the upper and lower bound deflections and slopes are improved. Bounds on the moment and shear agree with the exact result and are not plotted. Table 1 shows a comparison of the ten and fifty segment upper and lower bound transfer matrices over the entire length of the beam. Note that the bounds tend to converge with number of segments, but the number of computer operations is increasing. Table 2 gives the reader a brief comparison of the numerical results plotted in Figures 2 and 3 and the exact solution obtained from equation (8). Now consider the analysis of the statically indeterminate beam shown in Figure 4 under the loads P,, P2, . . . P,. The difference between this and the previous analysis is that one must determine bounds on the initial state vector before determining bounds on intermediate state vectors. The procedure for finding bounds on the initial state vector due to the q loads can be obtained by superposition of solution bounds for the various loads. For only P, applied to the structure, the boundary state vectors are related by
436
T. J. MCDANIEL AND V. R. MURTHY 0.E
0.5
0.4 z _; .i ; : ;_
0.:
0” 0.2
0-I
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
L-J
0.6
0.9
I.0
X/l
10 segments; ----,
Figure 2. Upper and lower bound deflections. -,
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.6
0.9
50 segments.
I.0
x/l
Figure 3. Upper and lower bound slopes. -,
10 segments; ----,
50 segments.
BOUNDS FOR VARYING GEOMETRY BEAMS
TABLE
431
1
Upper bounds (U.B.) and lower bounds (L.B.) for the elements of the transfer matrix compared with the exact result [T(l)] ; I = 6.5 in, EZ,, = 297.3132 lb-in2, CI= -1.0, P = 1 lb 10 Segments Element
f23
t24
L.B.
U.B.
0.1010633 0.2576360 0.0371916 0.1406822
0.1031049 0.2628406 0.0379429 0.1435241
0.0972474 0.2746312 0.0357190 0.1349262
0.1074750 O-2736748 0.0394756 0.1491165
f13
t14
50 Segments L.B.
U.B.
Exact 0.1020725 0*2602210 0.9375660 0.1421065
TABLE 2
(a) Upper and lower bound deflections
(in);
(b) upper and lower bound slopes (rad)
(4 50 Segments
10 Segments Station Xl1 0.0 0.5 1.0
Exact
U.B.
L.B.
U.B.
L.B.
0.0 0.1212 o-4509
::;049 0.3584
0.0 0.1141 O-41 25
0.0 0.1109 0.3940
0.0 0.1125 04032
U.B.
L.B.
Exact
(b) 10 Segments
Station Xl1
lJ,B.
50 Segments
L.B.
0.0 0.5
0.0
0.0
0.0
0.0
0.0
0.0732
0.0615
0.0684
O-0660
0.0672
l-0
0.1216
0.0830
0.1059
0.0982
0.1021
Sick view
Figure 4. Statically indeterminate beam with varying geometry.
438
T. J. MCDANIEL AND V. R. MURTHY
where [T] is the transfer matrix over the length of the beam and [S] is the transfer matrix from station q to station N. Once the boundary conditions are specified, e.g., for a clampedclamped beam, one can solve for the initial state vector. For the present example one obtains
P,lA, -f23s14
+
t13s24
where A = t13 tz4 - t,, t23. Since only bounds on [T] and [S] are presumed known, these are used to obtain bounds on MO and V, and hence on the initial state vector. Since all transfer matrix elements are non-negative for the present problem, an upper bound to MO would occur when tz4s14/d is as large as possible and t,,s,,lA is as small as possible for P, )O. It follows that for upper bounds
where d^ = 2,, f24 - i,,i,, ) 0 and A = ~,,i,, - f141z3 ) 0. The fact that d^ and 2 are positive can be verified from the numerical computation. Corresponding lower bounds on MO and V, are obtained for P, ) 0 by interchanging v and A in all terms of equation (13). If P, is negative the left-hand side of equation (13) receives an interchange of v and A for both bounds.
4. BOUNDS
FOR MULTISPAN
STRUCTURES
For a multispan beam with varying geometry, as shown in Figure 5, the development of bounds for the transfer matrix over the various segments of beam is described in previous
Figure 5. Multispan beam structure.
sections. Transmitting across the various supports of the multispan structure adds a new complexity to the bounding transfer matrix analysis since the transfer matrix across a spring support may contain a negative element. Such is the case if the support resists deflection of the beam. The negative element results from the fact that deflection and shear force do not have the same positive direction on the positive face of the beam: i.e., on the face whose normal points in the positive x direction. The transfer matrix across a point of support is given by
where 1 and r denote left and right of the support. The k, and k, are, respectively, the deflection and rotation spring constants at point Q as shown by Figure 5. If k, and k, are known, then the upper and lower bounds on the transfer matrix across the point Q are identical.
BOUNDS FOR VARYING GEOMETRY BEAMS
439
The procedure for multiplying upper and lower bound transfer matrices developed in equation (12a) is a special case for which the transfer matrix elements are non-negative. In the present case, where elements of the transfer matrix may be either positive or negative, a more general procedure for multiplying transfer matrices is required. For example, bounds on the elements of transfer matrix [C] can be obtained from bounds on the elements of transfer matrices [A] and [B] where
[Cl = MIPI.
(14)
Let the upper bounds on [A] and [B] be written as
[A] = [A+]+ [A-],
[B] = [B+] + [B-l,
[k] = [A+]+ [k-l,
[B] = [ri+] + [B-l.
and the lower bounds as
The positive and negative elements of the bounds are separated into two distinct matrices. It is not difficult to show that
[e+]= [A+][B+]+ [A-][B-1, [e-l = [A’][&] + [A-][b’],
[d-l = [a+][&] + [k-]@‘I, [C+]= [k’][B’] + [A-][B-1,
where [C] = [C’] + [C-l,
[Cl = [C’] + [C-l.
Computation of the transfer matrix for the complete structure proceeds as before. In the computation of the initial state vector one needs to develop rules for division of one bound by another. The necessary logic for programming this operation is as follows. Let the numerator be b, the denominator be c, and the quotient be d; then for an upper bound on d choose (z=&/c’
if
8=&/C
if
t;, E > 0,
a=h/C if a= b/t & > 0, c^< 0,
6, c^< o, if 6 < o, E > o.
Similarly, the rules for obtaining the lower bound to d can be established. Computation of intermediate state vectors proceeds as in equation (14), except that [B] is replaced by a column matrix, i.e., the initial state vector, and [C] is replaced by the intermediate state vector. 5. BOUNDS ON DYNAMIC
RESPONSE
When a beam structure undergoes harmonic motion in time, the amplitude of the state vector satisfies equation (9) except that ai4 is no longer zero: i.e., now ai4 = pw”t(x)b(x), where p is the mass per unit volume and w is the frequency of the harmonic motion. Note that the aij are still non-negative. The previously developed procedure, i.e., choosing the largest or smallest value of each a,, in the interval 0 < x < x0 and solving the corresponding constant coefficient differential equations, yields the required upper and lower bounds for a given segment of beam. Once bounds to the transfer matrix have been obtained for each segment, one can proceed via equation (14), or in the absence of supports via equation (12), to obtain upper and lower bounds for the transfer matrix over the entire structure. From these bounds one can obtain
440
T. J. MCDANIEL
AND V. R. MURTHY
bounds on the frequency determinant. As an example, consider a clamped-clamped without supports for which the frequency determinant is A(d) = f,,
fz4- i,, fz3,
beam
A(&) = i,, i14 - i,, fz3,
where iij, J(o), iij, and l(o) are upper and lower bounds to the unknown tij and d(o). The above result is more complicated if some of the tij‘s are negative, which can occur for a structure with spring supports. The procedure for obtaining bounds in this case, by applying equation (14), is clear. The zeros of d(G) and d(ui) provide upper and lower bounds on the natural frequencies of the system. It is clear that since A(h) 2 A(h), the first zero of A(d) is an upper bound to the first natural frequency, a lower bound to the second natural frequency, and so on. The reverse is true of zeros of the lower bound A(h). Bounds on the mode shapes of the system can be constructed once bounds on the natural frequencies are obtained. Let G and Lj be the bounds to natural frequency o: i.e., Q > o z &. For the above example of a clamped-clamped single span beam, one is free to choose at the natural frequency an amplitude such as V(0) = 1.0 and solve for the remaining unknown in the state vector M(0). In terms of bounds MO) = -~,,(~)/M~),
k(O) = -f14(oj)/i13(0).
From the above bounds on the state vector, one can find bounds on the mode shape. The upper deflection mode shape for the above case is given by a(x) = -i,,(x, ~2)~,,(h)/f,,(d) + i14(x,6). The lower bound is obtained by interchanging the upper and lower bound symbols in the above equation. The delta matrix method is an alternative approach to calculating the frequency determinant for structures analyzed by a transfer matrix approach. Basically, the delta matrix method reverses the order of computing the transfer matrix product and taking determinants. This tool is used to improve the accuracy of computation of natural frequencies. The rules for setting up the delta matrix method are given in reference [2] for the case where the transfer matrix elements are obtained from an analytical solution. The method is easily extended to the case where only bounds on the transfer matrix are known. For example, for a given segment of beam the one-one element of the delta matrix, i.e., AlI, is given by A,, = tll t,, t,, t,,. Bounds on Al, are given by
for a typical segment of beam. Similarly, bounds on the remaining elements of the delta matrix for the given segment can be obtained. Using the rules for obtaining bounds for matrix products given by equation (14), one can obtain bounds on the delta matrix for the entire structure, including spring supports, as a matrix product of delta matrix bounds for the various segments. Upper and lower bounds on the natural frequencies of a clamped-clamped single span beam with varying geometry were computed by the delta matrix procedure described above. As an example, bounds on the natural frequencies of the beam shown in Figure 1 are obtained. This case was chosen because the exact solution is available from Cranch and Alder [24] for comparison. Table 3 shows the comparison of the bounds with the exact result.
441
BOUNDSFORVARYING GEOMETRY BEAMS TABLE3 Bounds on natural frequencies;
I = 6-5 in, El0 = 808-l 811 lb+?,
p0 tobo = O-2005 x 10-l
lb-s2/in2, CI= -1.0 U.B. @ad/s)
Exact solution @ad/s)
LB. (rad/s) 300 Segments 600 Segments
Frequency
300 Segments
01
22.6956
22.5892
22.5004
22.2882
22.3858
-
142.9715
115+rO86
-
105.2926
02
600 Segments
As one can see from the above examples, the number of computer operations can become excessive as the accuracy of the bounds is increased. The efficient generation of transfer matrix bounds is essential and is given further consideration below. 6. IMPROVED BOUNDS ON THE TRANSFER MATRIX
As an alternative to the procedure of constructing transfer matrix bounds over a number of short segments of beam, the Picard iteration procedure is used to improve a coarse bound such as that given by equation (11) over the entire length of the beam. One can show that the governing differential equation for the beam, i.e., equation (lo), can be converted into the integral equation [T(x)] = j [.4(a)][T(i)] df + [I]. 0
The Picard iteration is a procedure for generating a sequence of approximations [T,(x)], [T,Wl, ... [T,(x)] which converge, provided [A(x)] is continuous, to [T(x)] as n approaches infinity. The [T,(x)]‘s are generated by (15)
Suppose one starts with [T,(x)] = [pl(x)], w h’ich is a coarse upper bound to the transfer matrix for the entire beam. Since the elements of [A(x)] are non-negative for static beams, it follows that the resulting [p2(x)] obtained from equation (15) is also an upper bound.Using differential inequalities one can show that [F2] is a better upper bound since it is also a lower bound to @‘dx)l.W’,(x)1 = [FCx)I IS ’ chosen to be a coarse lower bound, it follows by similar arguments that the [p2(x)] obtained by equation (15) is a better lower bound. The number of digits of agreement of the upper and lower bounds can be increased simply by increasing the number of Picard iterations which are carried out analytically. One finds
This iterative procedure reduces the number of computer operations since [T,(x)] can be computed analytically and finally evaluated numerically for use in determining bounds on the solution state vector. As an example, the upper bound transfer matrix for a beam with exponential width variation given by equation (11) was Picard iterated to obtain better upper bounds. After two iterations all elements of the upper bound [~Jx)] agreed with those of the exact solution [T(x)]. To illustrate the iterative technique a constant in-plane tension is introduced into the variable geometry beam. A closed form solution for the transfer matrix [T(x)] is not known
442
T. J. MCDANIEL AND V. R. MURTHY
in this case. The governing differential equation is now given by equation To obtain the first approximation to the upper and lower bound transfer [f3] to [T(x)] one solves
(O), with uJ2 = P. matrices [r,] and
and
where [d] minimum analytical iterations
> [A(x)] 3 [k]. The 2 is a constant coefficient matrix value in 0 < x < 1. k contains the maximum value of expressions for the elements of the bound transfer are given below. An exponential variation in bending
where El(x) is replaced by its El(x) in the given range. The matrix [7’3] after two Picard stiffness is considered :
El(x) = EIoexp(crxi/). For the first approximation,
El(x) is replaced by CEI,. After two Picard iterations
one obtains
where t,z = x + (Pl2Elo /NC, - C,) + PxIS1- Z~cr/ls,, t,3 = (1/2EIr~)(Cl + C,) +
CrX/lS,
-
S3/S3,
t,, = (1 /2EZ0 j?)(C, - C,) + x/S1 - 2ci/I&, 122 =
1 +
U/2-&
tz3 =
(1/2&)(c3
f32 =
Pf12,
I34 =
t12,
B)(C, + t33 =
-
c4)
C4) +
u/s,
+
P/S,, 1,
1 + (P/2EI,)(C,
t24 =
(1 WLJ
P)(C,
-
C4)
+
1 iSI,
+ C2) + PCrX/ZS, - PS3/S3
By choosing the value of C appropriately, one can obtain the upper and lower bound transfer matrices from the above expressions for the elements of the transfer matrices. For example, if a ) 0, then by choosing C = 1 the upper bound transfer matrix is obtained and by choosing C = exp(a) the lower bound transfer matrix is obtained for the entire length of the beam. The accuracy of upper and lower bound transfer matrices can be further improved by dividing the beam into N segments. For each segment the upper and lower bound transfer matrices can be obtained and improved by using the Picard iteration procedure. After obtaining the
t34
t24
t14
0.2736975 x 1O-3 0.1074890 x 1O-3 6.527370
/ 2 Iterations
TABLE 4
0.2735763 x 1O-3 0.1073930 x 10-a 6.527358
2 Iterations,
, 10 segments 0.2735760 x 1O-3 0.1073929 x 1O-.3 6.527358
200 steps
Runge-Kutta
T = 100.0 lb, I,+(~,= 0.01
10 segments
0.2735753 x 1O-3 0.1073927 x 1O-3 6.527358
I 2 Iterations,
transfer matrix eiements. I = 6.5 in, El0 = 1.041666 x 10’ lb-in2, cx= 1.0,
U.B. .A
Bounds on
L.B. A
\
0.2734032 x 1O-3 0.1072912 x 1O-3 6.527340
2 Iterations
U.B.
0.0 -0.3807096 -0.3441979 -0.2876234 04640264 0.1060432
0.0 0.2 0.4 0.6 0.8 1.0
,
0.0 -0.3016743 -0.8097403 -0.1076248 -0.8080679 +0.1748209
U.B.
xl1
Station
0.0 0.2 0.4 0.6 0.8 1.0
Station
x x x x x
x x x x x
lO-2 1o-2 1O-3 10-2 10-l
lo-’ IO-2 10-l 1O-2 1O-2
L.B.
o-0 -0.3974145 -0.3826083 -0.9237368 0.3727525 0.09397775
1 Segment
x x x x x
x x x x x
L.B.
0.0 -0.3118871 -0.8553294 -0.1187836 -1.020109 -0.1748209
1 Segment
1o-2 1O-2 1O-3 1O-2 10-l
1O-2 1O-2 10-l 1O-2 1O-2
\
’
and lower bound moments (in-lb);
c U.B.
0.0 -0.3891095 -0.3634552 -06061063 0.4183239 0.09998691
10 Segments h
0.0 -0.3890680 -03633622 -06046154 0.4185302 0~1000131
lO-2 1O-2 1O-3 1O-2 10-l
1O-2 1O-2 10-l lo-” 1O-5
L.B. x x x x x
x x x x x
L.B. 0.0 -03068206 -0.8326415 -0.1132210 -0.9143181 -04008712
10 Segments A
U.B.
,
0.0 -0.3067958 -0.8325307 -0.1131943 -0.9138206 +04008712
(4
x x x x x
x x x x x
1O-2 1O-3 1O-3 lO-2 10-l
l0-L 1O-2 10-l lO-2 1O-5
\
,
(c) upper (d) upper and lower bound shear force (lb)
TABLE 5 (a) Upper and lower bound dejections (in); (b) upper and lower bound slopes (rad);
R-K x x x x
0.0 -0.389088 x -0.3634087 x -0.6053601 x 0.4184261 x 0.1 x IO-’
R-K
0.0 -0.3068082 -0.8325860 -0.1132076 a.9140692 0.0
1O-2 1O-2 1O-3 lo-’
1O-2 1O-2 lo-’ 1O-2
U.B.
289.2014
All
I 283.1185
L.B.
-529.1155 -161.3734 206.1372 573.8588 942.0805 1310.980
-518.2715 -142.6113 232.8425 608.5379 984.7679 1361.713
1 Segment A
L.B.
1 Segment
U.B.
w/l
Station
0.0 0.2 0.4 0.6 0.8 1.0
Station
7
\
r
(4
I
286.2938
U.B.
286.2802
L.B.
223.7796 -151.9221 219.7164 591.5812 963.9634 1337.042
-523.7553 -151.8801 219.7761 591.6587 964.0589 1337.155
10 Segments *
L.B.
10 Segments * U.B.
TABLE 5-continued
,
\
286.2867
R-K
-523.7667 -151.9008 219.7462 591.6195 964.0102 1337.097
R-K
T. J. MCDANIEL
446
AND V. R. MURTHY
upper and lower bound transfer matrices for each segment, the transfer matrices for the whole system can be obtained by successive multiplications of transfer matrix bounds using equation (14). To illustrate the effect of Picard iteration, bounds on the state vector for a clamped-rotated clamp ($(I) given) beam with in-plane tension are computed. The variation of H(x) is again exponential. The present problem is analyzed by using two Picard iterations for the entire length of the beam and by dividing the beam into ten segments and again performing two Picard iterations on the solution. Selected elements of the upper and lower bound transfer matrices for the entire length of the beam are presented in Table 4. The upper and lower bound state vectors are presented in Table 5. Since the exact solution for this problem is not known, the results are compared with solutions obtained by the Runge-Kutta procedure. The difference between upper and lower bounds obtained in Table 5 cannot be distinguished on a plot. The actual time of computation on digital computers for obtaining upper and lower bound transfer matrices and state vectors at ten stations is significantly less than that required for ordinary transfer matrix analysis using Runge-Kutta. The comparison is shown below. Upper and lower bounds n \ 10 Segments
c 1 Segment 0.58 s
Transfer matrix analysis using Runge-Kutta
1.11 s
5.56 s
7. CONCLUSIONS In the above study a procedure for obtaining upper and lower bounds for linear boundary value problems is developed. As an application, beams with varying geometry are considered. These bounds are obtained by converting the boundary value problem into an initial value problem via a transfer matrix formulation. This approach allows the direct application of the theory of differential and integral inequalities, which is primarily applicable to the construction of bounds for initial value problems. A procedure for providing initial bounds on the transfer matrix and one for refining these bounds are presented. In particular, it is found that the Picard iteration is an effective procedure for refining the bounds. This procedure yields analytical bounds for the elements of a transfer matrix. These upper and lower bounds to elements of a transfer matrix agree to several places when evaluated numerically for a given length of the structure. It should be pointed out that even though exponential variation of El(x) was assumed for all examples considered, any piecewise continuous and bounded function could be used in its place. Many simple El(x) variations lead to quite complicated integrations, since it is the reciprocal of this function which must be integrated in the Picard iteration procedure. A means of avoiding this problem is to use a simple function such as a polynomial or an exponential function to provide simple upper and lower bounds to the function 1/El(x) in a given range. If these bounds are close, the bounds constructed on the simpler problem by the Picard iteration procedure will converge to close bounds on the actual problem. In the above analysis the upper and lower bounds on the transfer matrix elements for a static or dynamic beam are constructed independently. When elastic or rigid supports are included as a part of the structure, these bounds become coupled, Coupling of upper and lower bounds can occur for the beam alone if negative elements occur in the [A] matrix : e.g., for a beam-column analysis. From the transfer matrix bounds, one can proceed to the development of both upper and lower bounds to the entire state vector. Bounds on all
BOUNDSFOR VARYINGGEOMETRYBEAMS
447
quantities of the state vector can be obtained to essentially the same degree of precision. It is shown that the actual time of computation for generating the upper and lower bounds is significantly less than that required for an ordinary transfer matrix analysis using the RungeKutta procedure. The accuracy of the solution is apparent when upper and lower bounds are available, whereas this is not the case for the numerical solutions. It is necessary in the ordinary transfer matrix approach to use techniques such as the delta matrix method [2], formulation of a super matrix equation [9], etc., to circumvent numerical difficulties associated with the analysis of large and multispan structures. Such techniques can also be employed when bounds on the transfer matrix elements are available. It is important to use all possible information that the differential inequalities can supply in order that the above techniques will give bounds which are sufficiently accurate to be of use to engineers. ACKNOWLEDGMENTS This research was sponsored by the National Science Foundation under Grant No. GK-40589 and by the Iowa State University Engineering Research Institute. Thanks are due to Mr P. V. Thangam
Babu and Miss Cecelia Byers for assistance
with this manuscript.
REFERENCES 1. T. J. MCDANIEL and J. P. HENDERSON1974 The Shock and Vibrution Digest 6-l (Part I), -2 (Part II). Review of transfer matrix vibration analysis of skin stringer structure. 2. E. C. PESTELand F. A. LECKIE1963 Matrix Methods in Elastomechanics. New York: McGrawHill Book Company, Inc. 3. Y. K. LIN 1967 Probabilistic Theory of Structural Dynamics. New York: McGraw-Hill Book Company, Inc. 4. Y. K. LIN and T. J. MCDANIEL1969 Journal of Engineering for Industry 91,1133-l 141. Dynamics of beam-type periodic structures. 5. C. A. MERCERand C. SEAVY1967 Journal of Sound and Vibration 6,149-l 62. Prediction of natural frequencies and normal modes of skin-stringer panel rows. 6. D. J. MEAD 1971 Journal of Engineering for Industry 93, 783-792. Vibration response and wave propagation in periodic structures. 7. D. J. MEAD and G. SENGUPTA 1970 USAF Report AFML-TR-70-13. Propagation of flexural waves in infinite, damped rib-skin structures. 8. J. P. HENDERSON1972 Ph.D. Thesis, Ohio State University. Vibration analysis of curved skinstringer structures having tuned elastomeric dampers. 9. J. P. HENDERSONand T. J. MCDANIEL 1971 Journal of Sound and Vibration 18, 203-219. The analysis of curved multi-span structures. 10. T. J. MCDANIEL 1971 Journal of Aircraft 8, 143-149. Dynamics of circular periodic structures. 11. T. J. MCDANIEL 1792 Journal of Sound and Vibration 23, 217-227. Dynamics of non-circular stiffened cylindrical shells. 12. W. F. HUNTER 1970 NASA TN-D 6064. Integrating matrix method for determining the natural vibration characteristics of propeller blades. 13 T. J. MCDANIEL and J. D. LOGAN 1971 Journal of Sound and Vibration 19, 39-48. Dynamics of cylindrical shells with variable curvature. 14. W. WALTER1970 Diflerential and Integral Inequalities. New York: Springer Verlag. 15. M. H. PaorrEa and H. S. WEINBERGER1967 Maximum Principles in D@rential Equcttions. New Jersey: Prentice-Hall, Inc. 16. V. LAKSHMIKANTHAM and S. LEELA1969 D@erentialandIntegral Inequalities, Volume I, Ordinary Differential Equations. New York and London: Academic Press. 17. B. FRAEJIS DE VEUBEKE 1965 Proceedings of Conference on Matrix Methods in Structural Mechanics, AFFDL-TR-66-80, 863-886. Bending, and stretching of plates-special models for upper and lower bounds. 18. N. W. BAZLEY, D. W. Fox and J. T. STADTER1967 Zeitschrift fur angewandte Mathematik und Mechanik 47, 191-198. Upper and lower bounds for the frequencies of rectangular clamped plates.
448
T. I. MCDANIELAND V. R. MURTHY
19. K. NICKEL 1962 Ingenieur Archiv 31, 85-l 00. Eine einfache Abschatzung fur Grenzschichten. 20. E. ADAMS 1967 Acta Mechanica 3, 268-277. Rigorous bounds of exact solutions of boundary layer problems of Prandtl type. 21. F. JOHN 1965 Courant Institute of Mathematical Sciences, New York University, Research Report No. IMM-NYU 308. Estimates of the error in the equations of non-linear plate theory. 22. J. R. CANNON 1966 Proceedings of an Advanced Symposium Conducted by the Mathematics Research Centre, United States Army, at the University of Wisconsin, Madison. 21-54. Numerical solutions of nonlinear differential equations. 23. C. A. BELL and F. C. APPL 1974 JournalofAppliedMechanics 40, 1097-l 102. Bounds for initial value problems. 24. E. T. CRANCH and A. A. ALDER 1955 Journal of Applied Mechanics 21, 103-108. Bending vibrations of variab!e section beams.
SOLUTION
BOUNDS
FOR
T. J. MCDANIEL
VARYING AND
GEOMETRY
BEAMS
V. R. MURTHY
Department of Aerospace Engineering and Engineering Research Institute, Iowa State University, Ames, Iowa 50010, U.S.A. (Received 7 April 1975, and in revised form 15 July 1975)
The theory of differential and integral inequalities is applied to obtain upper and lower bounds to the transfer matrix for beams with varying geometry. Various techniques of generating and refining these bounds are investigated. Numerical results indicate that these bounds can be refined to produce numerical agreement of the upper and the lower bound to a given number of significant digits. Proceeding from bounds on the transfer matrix elements a theory is developed for determining upper and lower bounds on the natural frequencies and mode shapes and on the solution state vector for static loading of such beams. This procedure is then extended to the analysis of multispan beams with varying geometry. Numerical results are presented for various configurations.
1. INTRODUCTION
A recent survey article by McDaniel and Henderson [l] indicates the variety of multispan structures that have been analyzed by the transfer matrix method. The theoretical and numerical aspects of transfer matrix analysis have been developed by a number of researchers, including Pestel and Leckie [2], Lin [3, 41, Mercer [5], Mead [6, 71, Henderson [8, 93, and McDaniel [lo, 111, for application to the dynamics of stiffened beams, plates, and shells. Structures with varying geometry, i.e., varying thickness, varying curvature, etc., can be considered in the transfer matrix analysis by approximating the varying geometry structure by constant geometry segments or by using numerical techniques such as Runge-Kutta, predictor-corrector, or the integrating matrix procedure [12] to obtain numerical transfer matrices. Due to the difficulty of solving variable coefficient differential equations associated with the varying geometry structure, few exact solutions are available for comparison with these techniques. An exact solution is known for simply supported cylindrical shells with varying curvature only for the case of exponentially varying curvature [ 131. The use of small integration steps and double precision arithmetic is standard practice for generating numerical transfer matrices. Even so, it is difficult to determine the accuracy of the numerical method and to detect cases where the accuracy of the numerical method is insufficient. It is important that the numerical procedures have sufficient accuracy but not excessive accuracy, as this increases significantly the cost of computation. If limited accuracy computations are used to determine the individual transfer matrices for a single or multispan structure, how does this affect the final solution? In the following sections this problem will be explored by developing upper and lower bounds for the individual transfer matrices. These bounds will then be used to determine bounds on the final solution and thus to determine the effect of limited accuracy intermediate computations on the final results. 431
T. J. MCDANIEL
432
AND V. R. MURTHY
2. DIFFERENTIAL
EQUALITIES
The theory of differential inequalities can be used to establish bounds on the solution of differential and integral equations and thus can provide answers to the above questions. The texts of Walter [14], Protter and Weinberger [I 51, and Lakshmikantham and Leela [ 161 summarize the main results and provide numerous references. For engineering application, the capability of establishing both upper and lower bounds through the theory of inequalities is quite important. It allows one to bracket unknown exact solutions by rough upper and lower bounds if this is sufficient and by refined upper and lower bounds which agree to a prescribed number of digits when required. Unlike energy methods which insure that complete continuity of the displacement field will produce lower bounds on the static influence coefficients [17] or that complete continuity of the stress field in a structure provides upper bounds to the static influence coefficients, the method of differential and integral inequalities provides both upper and lower bounds to the transfer matrix elements for both the static or harmonically vibrating structure. This suggests that it may be possible to construct upper and lower bounds on natural frequencies and on the corresponding mode shapes whereas the Rayleigh-Ritz and the method of Bazley, Fox, and Stadter [18] yield bounds on only natural frequencies. The primary engineering application of the theory of differential and integral inequalities is to construct bounds on linear and non-linear initial value problems. Since the transfer matrix approach to structural analysis converts the two point linear boundary value problem into an initial value problem, the inequality theorems apply directly. The inequality theorems have not been developed sufficiently to allow the construction of bounds for non-linear, two point boundary value problems except for special classes of problems. Even though extensive development of the theoretical aspects of differential inequalities can be found, few engineering application papers are accompanied by such bounds. Some applications, however, can be cited such as those of Nickel [19] and Adams [20], who have explored boundary layer flow equations using this approach. John [21] has developed error estimates for a non-linear plate analysis, and Cannon [22] has found bounds to a heat conduction problem. In general, high precision in the computed results was not demonstrated. However, in a recent article Bell and Appl [23] developed a computational procedure based on differential inequalities which provided quite accurate bounds. Upper and lower bounds for the response of a non-linear single degree of freedom system agreed to several digits over the range of computation for cases where the system was excited by short duration excitation or provided with given initial conditions. The differential inequalities used in the present beam analysis are of the following form. Consider a linear initial value problem :
;
u-(x)1 = [4-~)lmx)l,
O
(1)
with initial conditions
P-m = [II, and with the elements
of [A(x)], i.e., aij(x), satisfying L7ij
If it is assumed
that aij is piecewise continuous
>
0.
and bounded,
(3) then existence
and uniqueness
BOUNDS FOR VARYING GEOMETRY BEAMS
433
of the solution is assured. Then if [F(X)] and [i(x)] satisfy the following differential inequalities,
wmx)l,
m’(o)12 VI,
(4)
G Lw1m1,
[WV G VI,
(5)
& VP)1 2
ml
$ it follows that
ml
2
.,
Lwl
2
mx)l.
(6)
A proof of this result can be found in Walter’s book [14]. In fact, the theorems presented in this text apply to the case where the elements of [A] may be negative and to non-linear differential equations of the initial value type. Since the elements of [A] are non-negative in the present static or dynamic beam analysis, it can be shown that the solution [T(x)] of equations (1) and (2) has non-negative elements. This can be shown by representing the solution [T(x)] in terms of the matrizant of [A(x)]. This fact is quite useful in the construction of the inequalities (4) and (5). 3. BOUNDS FOR THE STATIC DEFLECTION OF BEAMS One can apply the differential inequalities discussed in the last section to determine bounds on the behavior of beams with variable geometry. Since an exact solution for the transfer matrix for a beam with exponential width variation and constant thickness shown in Figure 1
Figure 1. Beam with exponentialvariation of width. b(x) = b,,exp(ax/l); t(x) = rO; E(x) = Eo. easily obtained, this solution can be compared, as a first example, with bounds on the transfer matrix elements. The exact solution for the transfer matrix is given by
is
‘1
x
t,,
t,, -
1
l23
?24
0
001
x
000
1
(7)
where -ax/l}/a - l/a) I/EZc f&) = (x + t14(x) = (1 + exp {-cl.x/1})Z2x/EZo ~1~+ (exp {-ax/f} - 1) 2f3/EZoa3, [23(x)
=
(1
t24(x)
=
(l/a
-
exp -
{-a//})
(l/a
El0 = Eobo t,3,‘12.
+
&f.?z,, x)
exp
tl,
{-Rx//})
&??&
a,
434
T. J. MCDANIEL
AND V. R. MURTHY
For the cantilever beam with a tip load P shown in Figure analysis that the state vector at station x is given by
WCC, II/(Xl
1, one finds via the transfer
.t,31-t,,
III 1 M,,,
=
matrix
t2J - t24 l-x
V(X)
P.
(8)
-1
where W, $I,Mand Vare usual notation for the deflection, slope, moment, and shear contained in the state vector for a beam. Equation (8) is obtained by substituting boundary conditions W(0) =$(O) = M(I) = 0 and V(I) = -P in the transfer matrix equation (7), solving for the state vector at x = 0, and transmitting that state vector to station x. The initial state vector can be obtained in this particular case from statics and is precisely known. Now, let us apply differential inequalities to obtain bounds on the state vector as a function of x. A comparison of bounds can be made with the exact solution only in special examples such as the present one. In most cases one will be able to compare only the upper and lower bounds. To develop bounds on the transfer matrix, one proceeds from the governing differential equation
Recall that the transfer
with the initial which satisfies
condition
matrix also satisfies the same diherential
equation,
[T(O)] = [I]. Let us choose an approximate
& [Rx)] 2 [A(x)][Qx>]
in
i.e.,
transfer
matrix
[T(x)]
o
and [p(O)] = [r(O)] = [I]. A ccording to the previous section it follows that [F(X)] is an upper bound to [T(x)] for the specified range of x. It is easily shown for the case u ) 0 that the solution of equation (10) with CL= 0 yields a [p(x)] which satisfies the above conditions and is given by t^13= x2/2EI,,
(11) iz4 = x2/2EI,. This upper bound transfer matrix is clearly generated by choosing the smallest El(x) in the given interval. It follows that a lower bound transfer matrix [T(x)] is obtained by choosing the largest value of El(x), i.e., EI,exp{crx,/l}, in the given interval and solving the resulting
BOUNDS
FOR VARYING
GEOMETRY
435
BEAMS
constant coefficient differential equation. This lower bound satisfies the required conditions
and [f(O)] = [T(O)] = [I]. The advantage of the present method of choosing the upper and lower bound transfer matrices is that the required inequalities are easily shown to be satisfied. The disadvantage is that if El(x) varies rapidly along the structure, the range over which [F] and [f] are acceptable bounds is considerably reduced: i.e., x0 is small. A coarse bound on {Z(x)} for the present example is obtained by letting x,, = 1 so that [fWl and PI” x )I are bounds over the entire length of the beam. The resulting boundary state vectors (p(x)} and {k(x)} are given by 1x2/2EZ0- x3j6EZo for
u > 0.
(12)
{i(x)} is obtained in this case by replacing El0 by EZoexp{ctx,/l} in equation (12). One can easily verify that {z} and {i} are upper and lower bounds to {Z} in equation (8). If c1is negative then the upper and lower bounds of equation (12) switch roles. The bounds obtained are coarse since the transfer matrix bounds from which they are derived are not particularly close bounds. Improved upper and lower bounds can be constructed by subdividing the beam into N segments and finding upper and lower bounds to the transfer matrix in each interval. Since for the static or dynamic analysis of beams the elements of a transfer mtarix are non-negative, i.e., tij 2 0, it follows that the upper and lower transfer matrices overj segments are given by PI = rf+,lv+~-Il~~~ rf+21m,
[?I = [f,][i;_,].
. . [7’23[TJ.
(124
As in the previous example, the maximum and minimum values of q in each interval were used to construct the upper and lower bound transfer matrices for that interval. The beam of Figure 1 was divided into ten segments, and the resulting bounds on deflection and slope are plotted as solid lines in Figures 2 and 3. Dividing the beam into fifty segments produced the bounds on deflection and slope given by dashed lines in the same figures. In the latter case the upper and lower bound deflections and slopes are improved. Bounds on the moment and shear agree with the exact result and are not plotted. Table 1 shows a comparison of the ten and fifty segment upper and lower bound transfer matrices over the entire length of the beam. Note that the bounds tend to converge with number of segments, but the number of computer operations is increasing. Table 2 gives the reader a brief comparison of the numerical results plotted in Figures 2 and 3 and the exact solution obtained from equation (8). Now consider the analysis of the statically indeterminate beam shown in Figure 4 under the loads P,, P2, . . . P,. The difference between this and the previous analysis is that one must determine bounds on the initial state vector before determining bounds on intermediate state vectors. The procedure for finding bounds on the initial state vector due to the q loads can be obtained by superposition of solution bounds for the various loads. For only P, applied to the structure, the boundary state vectors are related by
436
T. J. MCDANIEL AND V. R. MURTHY 0.E
0.5
0.4 z _; .i ; : ;_
0.:
0” 0.2
0-I
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
L-J
0.6
0.9
I.0
X/l
10 segments; ----,
Figure 2. Upper and lower bound deflections. -,
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.6
0.9
50 segments.
I.0
x/l
Figure 3. Upper and lower bound slopes. -,
10 segments; ----,
50 segments.
BOUNDS FOR VARYING GEOMETRY BEAMS
TABLE
431
1
Upper bounds (U.B.) and lower bounds (L.B.) for the elements of the transfer matrix compared with the exact result [T(l)] ; I = 6.5 in, EZ,, = 297.3132 lb-in2, CI= -1.0, P = 1 lb 10 Segments Element
f23
t24
L.B.
U.B.
0.1010633 0.2576360 0.0371916 0.1406822
0.1031049 0.2628406 0.0379429 0.1435241
0.0972474 0.2746312 0.0357190 0.1349262
0.1074750 O-2736748 0.0394756 0.1491165
f13
t14
50 Segments L.B.
U.B.
Exact 0.1020725 0*2602210 0.9375660 0.1421065
TABLE 2
(a) Upper and lower bound deflections
(in);
(b) upper and lower bound slopes (rad)
(4 50 Segments
10 Segments Station Xl1 0.0 0.5 1.0
Exact
U.B.
L.B.
U.B.
L.B.
0.0 0.1212 o-4509
::;049 0.3584
0.0 0.1141 O-41 25
0.0 0.1109 0.3940
0.0 0.1125 04032
U.B.
L.B.
Exact
(b) 10 Segments
Station Xl1
lJ,B.
50 Segments
L.B.
0.0 0.5
0.0
0.0
0.0
0.0
0.0
0.0732
0.0615
0.0684
O-0660
0.0672
l-0
0.1216
0.0830
0.1059
0.0982
0.1021
Sick view
Figure 4. Statically indeterminate beam with varying geometry.
438
T. J. MCDANIEL AND V. R. MURTHY
where [T] is the transfer matrix over the length of the beam and [S] is the transfer matrix from station q to station N. Once the boundary conditions are specified, e.g., for a clampedclamped beam, one can solve for the initial state vector. For the present example one obtains
P,lA, -f23s14
+
t13s24
where A = t13 tz4 - t,, t23. Since only bounds on [T] and [S] are presumed known, these are used to obtain bounds on MO and V, and hence on the initial state vector. Since all transfer matrix elements are non-negative for the present problem, an upper bound to MO would occur when tz4s14/d is as large as possible and t,,s,,lA is as small as possible for P, )O. It follows that for upper bounds
where d^ = 2,, f24 - i,,i,, ) 0 and A = ~,,i,, - f141z3 ) 0. The fact that d^ and 2 are positive can be verified from the numerical computation. Corresponding lower bounds on MO and V, are obtained for P, ) 0 by interchanging v and A in all terms of equation (13). If P, is negative the left-hand side of equation (13) receives an interchange of v and A for both bounds.
4. BOUNDS
FOR MULTISPAN
STRUCTURES
For a multispan beam with varying geometry, as shown in Figure 5, the development of bounds for the transfer matrix over the various segments of beam is described in previous
Figure 5. Multispan beam structure.
sections. Transmitting across the various supports of the multispan structure adds a new complexity to the bounding transfer matrix analysis since the transfer matrix across a spring support may contain a negative element. Such is the case if the support resists deflection of the beam. The negative element results from the fact that deflection and shear force do not have the same positive direction on the positive face of the beam: i.e., on the face whose normal points in the positive x direction. The transfer matrix across a point of support is given by
where 1 and r denote left and right of the support. The k, and k, are, respectively, the deflection and rotation spring constants at point Q as shown by Figure 5. If k, and k, are known, then the upper and lower bounds on the transfer matrix across the point Q are identical.
BOUNDS FOR VARYING GEOMETRY BEAMS
439
The procedure for multiplying upper and lower bound transfer matrices developed in equation (12a) is a special case for which the transfer matrix elements are non-negative. In the present case, where elements of the transfer matrix may be either positive or negative, a more general procedure for multiplying transfer matrices is required. For example, bounds on the elements of transfer matrix [C] can be obtained from bounds on the elements of transfer matrices [A] and [B] where
[Cl = MIPI.
(14)
Let the upper bounds on [A] and [B] be written as
[A] = [A+]+ [A-],
[B] = [B+] + [B-l,
[k] = [A+]+ [k-l,
[B] = [ri+] + [B-l.
and the lower bounds as
The positive and negative elements of the bounds are separated into two distinct matrices. It is not difficult to show that
[e+]= [A+][B+]+ [A-][B-1, [e-l = [A’][&] + [A-][b’],
[d-l = [a+][&] + [k-]@‘I, [C+]= [k’][B’] + [A-][B-1,
where [C] = [C’] + [C-l,
[Cl = [C’] + [C-l.
Computation of the transfer matrix for the complete structure proceeds as before. In the computation of the initial state vector one needs to develop rules for division of one bound by another. The necessary logic for programming this operation is as follows. Let the numerator be b, the denominator be c, and the quotient be d; then for an upper bound on d choose (z=&/c’
if
8=&/C
if
t;, E > 0,
a=h/C if a= b/t & > 0, c^< 0,
6, c^< o, if 6 < o, E > o.
Similarly, the rules for obtaining the lower bound to d can be established. Computation of intermediate state vectors proceeds as in equation (14), except that [B] is replaced by a column matrix, i.e., the initial state vector, and [C] is replaced by the intermediate state vector. 5. BOUNDS ON DYNAMIC
RESPONSE
When a beam structure undergoes harmonic motion in time, the amplitude of the state vector satisfies equation (9) except that ai4 is no longer zero: i.e., now ai4 = pw”t(x)b(x), where p is the mass per unit volume and w is the frequency of the harmonic motion. Note that the aij are still non-negative. The previously developed procedure, i.e., choosing the largest or smallest value of each a,, in the interval 0 < x < x0 and solving the corresponding constant coefficient differential equations, yields the required upper and lower bounds for a given segment of beam. Once bounds to the transfer matrix have been obtained for each segment, one can proceed via equation (14), or in the absence of supports via equation (12), to obtain upper and lower bounds for the transfer matrix over the entire structure. From these bounds one can obtain
440
T. J. MCDANIEL
AND V. R. MURTHY
bounds on the frequency determinant. As an example, consider a clamped-clamped without supports for which the frequency determinant is A(d) = f,,
fz4- i,, fz3,
beam
A(&) = i,, i14 - i,, fz3,
where iij, J(o), iij, and l(o) are upper and lower bounds to the unknown tij and d(o). The above result is more complicated if some of the tij‘s are negative, which can occur for a structure with spring supports. The procedure for obtaining bounds in this case, by applying equation (14), is clear. The zeros of d(G) and d(ui) provide upper and lower bounds on the natural frequencies of the system. It is clear that since A(h) 2 A(h), the first zero of A(d) is an upper bound to the first natural frequency, a lower bound to the second natural frequency, and so on. The reverse is true of zeros of the lower bound A(h). Bounds on the mode shapes of the system can be constructed once bounds on the natural frequencies are obtained. Let G and Lj be the bounds to natural frequency o: i.e., Q > o z &. For the above example of a clamped-clamped single span beam, one is free to choose at the natural frequency an amplitude such as V(0) = 1.0 and solve for the remaining unknown in the state vector M(0). In terms of bounds MO) = -~,,(~)/M~),
k(O) = -f14(oj)/i13(0).
From the above bounds on the state vector, one can find bounds on the mode shape. The upper deflection mode shape for the above case is given by a(x) = -i,,(x, ~2)~,,(h)/f,,(d) + i14(x,6). The lower bound is obtained by interchanging the upper and lower bound symbols in the above equation. The delta matrix method is an alternative approach to calculating the frequency determinant for structures analyzed by a transfer matrix approach. Basically, the delta matrix method reverses the order of computing the transfer matrix product and taking determinants. This tool is used to improve the accuracy of computation of natural frequencies. The rules for setting up the delta matrix method are given in reference [2] for the case where the transfer matrix elements are obtained from an analytical solution. The method is easily extended to the case where only bounds on the transfer matrix are known. For example, for a given segment of beam the one-one element of the delta matrix, i.e., AlI, is given by A,, = tll t,, t,, t,,. Bounds on Al, are given by
for a typical segment of beam. Similarly, bounds on the remaining elements of the delta matrix for the given segment can be obtained. Using the rules for obtaining bounds for matrix products given by equation (14), one can obtain bounds on the delta matrix for the entire structure, including spring supports, as a matrix product of delta matrix bounds for the various segments. Upper and lower bounds on the natural frequencies of a clamped-clamped single span beam with varying geometry were computed by the delta matrix procedure described above. As an example, bounds on the natural frequencies of the beam shown in Figure 1 are obtained. This case was chosen because the exact solution is available from Cranch and Alder [24] for comparison. Table 3 shows the comparison of the bounds with the exact result.
441
BOUNDSFORVARYING GEOMETRY BEAMS TABLE3 Bounds on natural frequencies;
I = 6-5 in, El0 = 808-l 811 lb+?,
p0 tobo = O-2005 x 10-l
lb-s2/in2, CI= -1.0 U.B. @ad/s)
Exact solution @ad/s)
LB. (rad/s) 300 Segments 600 Segments
Frequency
300 Segments
01
22.6956
22.5892
22.5004
22.2882
22.3858
-
142.9715
115+rO86
-
105.2926
02
600 Segments
As one can see from the above examples, the number of computer operations can become excessive as the accuracy of the bounds is increased. The efficient generation of transfer matrix bounds is essential and is given further consideration below. 6. IMPROVED BOUNDS ON THE TRANSFER MATRIX
As an alternative to the procedure of constructing transfer matrix bounds over a number of short segments of beam, the Picard iteration procedure is used to improve a coarse bound such as that given by equation (11) over the entire length of the beam. One can show that the governing differential equation for the beam, i.e., equation (lo), can be converted into the integral equation [T(x)] = j [.4(a)][T(i)] df + [I]. 0
The Picard iteration is a procedure for generating a sequence of approximations [T,(x)], [T,Wl, ... [T,(x)] which converge, provided [A(x)] is continuous, to [T(x)] as n approaches infinity. The [T,(x)]‘s are generated by (15)
Suppose one starts with [T,(x)] = [pl(x)], w h’ich is a coarse upper bound to the transfer matrix for the entire beam. Since the elements of [A(x)] are non-negative for static beams, it follows that the resulting [p2(x)] obtained from equation (15) is also an upper bound.Using differential inequalities one can show that [F2] is a better upper bound since it is also a lower bound to @‘dx)l.W’,(x)1 = [FCx)I IS ’ chosen to be a coarse lower bound, it follows by similar arguments that the [p2(x)] obtained by equation (15) is a better lower bound. The number of digits of agreement of the upper and lower bounds can be increased simply by increasing the number of Picard iterations which are carried out analytically. One finds
This iterative procedure reduces the number of computer operations since [T,(x)] can be computed analytically and finally evaluated numerically for use in determining bounds on the solution state vector. As an example, the upper bound transfer matrix for a beam with exponential width variation given by equation (11) was Picard iterated to obtain better upper bounds. After two iterations all elements of the upper bound [~Jx)] agreed with those of the exact solution [T(x)]. To illustrate the iterative technique a constant in-plane tension is introduced into the variable geometry beam. A closed form solution for the transfer matrix [T(x)] is not known
442
T. J. MCDANIEL AND V. R. MURTHY
in this case. The governing differential equation is now given by equation To obtain the first approximation to the upper and lower bound transfer [f3] to [T(x)] one solves
(O), with uJ2 = P. matrices [r,] and
and
where [d] minimum analytical iterations
> [A(x)] 3 [k]. The 2 is a constant coefficient matrix value in 0 < x < 1. k contains the maximum value of expressions for the elements of the bound transfer are given below. An exponential variation in bending
where El(x) is replaced by its El(x) in the given range. The matrix [7’3] after two Picard stiffness is considered :
El(x) = EIoexp(crxi/). For the first approximation,
El(x) is replaced by CEI,. After two Picard iterations
one obtains
where t,z = x + (Pl2Elo /NC, - C,) + PxIS1- Z~cr/ls,, t,3 = (1/2EIr~)(Cl + C,) +
CrX/lS,
-
S3/S3,
t,, = (1 /2EZ0 j?)(C, - C,) + x/S1 - 2ci/I&, 122 =
1 +
U/2-&
tz3 =
(1/2&)(c3
f32 =
Pf12,
I34 =
t12,
B)(C, + t33 =
-
c4)
C4) +
u/s,
+
P/S,, 1,
1 + (P/2EI,)(C,
t24 =
(1 WLJ
P)(C,
-
C4)
+
1 iSI,
+ C2) + PCrX/ZS, - PS3/S3
By choosing the value of C appropriately, one can obtain the upper and lower bound transfer matrices from the above expressions for the elements of the transfer matrices. For example, if a ) 0, then by choosing C = 1 the upper bound transfer matrix is obtained and by choosing C = exp(a) the lower bound transfer matrix is obtained for the entire length of the beam. The accuracy of upper and lower bound transfer matrices can be further improved by dividing the beam into N segments. For each segment the upper and lower bound transfer matrices can be obtained and improved by using the Picard iteration procedure. After obtaining the
t34
t24
t14
0.2736975 x 1O-3 0.1074890 x 1O-3 6.527370
/ 2 Iterations
TABLE 4
0.2735763 x 1O-3 0.1073930 x 10-a 6.527358
2 Iterations,
, 10 segments 0.2735760 x 1O-3 0.1073929 x 1O-.3 6.527358
200 steps
Runge-Kutta
T = 100.0 lb, I,+(~,= 0.01
10 segments
0.2735753 x 1O-3 0.1073927 x 1O-3 6.527358
I 2 Iterations,
transfer matrix eiements. I = 6.5 in, El0 = 1.041666 x 10’ lb-in2, cx= 1.0,
U.B. .A
Bounds on
L.B. A
\
0.2734032 x 1O-3 0.1072912 x 1O-3 6.527340
2 Iterations
U.B.
0.0 -0.3807096 -0.3441979 -0.2876234 04640264 0.1060432
0.0 0.2 0.4 0.6 0.8 1.0
,
0.0 -0.3016743 -0.8097403 -0.1076248 -0.8080679 +0.1748209
U.B.
xl1
Station
0.0 0.2 0.4 0.6 0.8 1.0
Station
x x x x x
x x x x x
lO-2 1o-2 1O-3 10-2 10-l
lo-’ IO-2 10-l 1O-2 1O-2
L.B.
o-0 -0.3974145 -0.3826083 -0.9237368 0.3727525 0.09397775
1 Segment
x x x x x
x x x x x
L.B.
0.0 -0.3118871 -0.8553294 -0.1187836 -1.020109 -0.1748209
1 Segment
1o-2 1O-2 1O-3 1O-2 10-l
1O-2 1O-2 10-l 1O-2 1O-2
\
’
and lower bound moments (in-lb);
c U.B.
0.0 -0.3891095 -0.3634552 -06061063 0.4183239 0.09998691
10 Segments h
0.0 -0.3890680 -03633622 -06046154 0.4185302 0~1000131
lO-2 1O-2 1O-3 1O-2 10-l
1O-2 1O-2 10-l lo-” 1O-5
L.B. x x x x x
x x x x x
L.B. 0.0 -03068206 -0.8326415 -0.1132210 -0.9143181 -04008712
10 Segments A
U.B.
,
0.0 -0.3067958 -0.8325307 -0.1131943 -0.9138206 +04008712
(4
x x x x x
x x x x x
1O-2 1O-3 1O-3 lO-2 10-l
l0-L 1O-2 10-l lO-2 1O-5
\
,
(c) upper (d) upper and lower bound shear force (lb)
TABLE 5 (a) Upper and lower bound dejections (in); (b) upper and lower bound slopes (rad);
R-K x x x x
0.0 -0.389088 x -0.3634087 x -0.6053601 x 0.4184261 x 0.1 x IO-’
R-K
0.0 -0.3068082 -0.8325860 -0.1132076 a.9140692 0.0
1O-2 1O-2 1O-3 lo-’
1O-2 1O-2 lo-’ 1O-2
U.B.
289.2014
All
I 283.1185
L.B.
-529.1155 -161.3734 206.1372 573.8588 942.0805 1310.980
-518.2715 -142.6113 232.8425 608.5379 984.7679 1361.713
1 Segment A
L.B.
1 Segment
U.B.
w/l
Station
0.0 0.2 0.4 0.6 0.8 1.0
Station
7
\
r
(4
I
286.2938
U.B.
286.2802
L.B.
223.7796 -151.9221 219.7164 591.5812 963.9634 1337.042
-523.7553 -151.8801 219.7761 591.6587 964.0589 1337.155
10 Segments *
L.B.
10 Segments * U.B.
TABLE 5-continued
,
\
286.2867
R-K
-523.7667 -151.9008 219.7462 591.6195 964.0102 1337.097
R-K
T. J. MCDANIEL
446
AND V. R. MURTHY
upper and lower bound transfer matrices for each segment, the transfer matrices for the whole system can be obtained by successive multiplications of transfer matrix bounds using equation (14). To illustrate the effect of Picard iteration, bounds on the state vector for a clamped-rotated clamp ($(I) given) beam with in-plane tension are computed. The variation of H(x) is again exponential. The present problem is analyzed by using two Picard iterations for the entire length of the beam and by dividing the beam into ten segments and again performing two Picard iterations on the solution. Selected elements of the upper and lower bound transfer matrices for the entire length of the beam are presented in Table 4. The upper and lower bound state vectors are presented in Table 5. Since the exact solution for this problem is not known, the results are compared with solutions obtained by the Runge-Kutta procedure. The difference between upper and lower bounds obtained in Table 5 cannot be distinguished on a plot. The actual time of computation on digital computers for obtaining upper and lower bound transfer matrices and state vectors at ten stations is significantly less than that required for ordinary transfer matrix analysis using Runge-Kutta. The comparison is shown below. Upper and lower bounds n \ 10 Segments
c 1 Segment 0.58 s
Transfer matrix analysis using Runge-Kutta
1.11 s
5.56 s
7. CONCLUSIONS In the above study a procedure for obtaining upper and lower bounds for linear boundary value problems is developed. As an application, beams with varying geometry are considered. These bounds are obtained by converting the boundary value problem into an initial value problem via a transfer matrix formulation. This approach allows the direct application of the theory of differential and integral inequalities, which is primarily applicable to the construction of bounds for initial value problems. A procedure for providing initial bounds on the transfer matrix and one for refining these bounds are presented. In particular, it is found that the Picard iteration is an effective procedure for refining the bounds. This procedure yields analytical bounds for the elements of a transfer matrix. These upper and lower bounds to elements of a transfer matrix agree to several places when evaluated numerically for a given length of the structure. It should be pointed out that even though exponential variation of El(x) was assumed for all examples considered, any piecewise continuous and bounded function could be used in its place. Many simple El(x) variations lead to quite complicated integrations, since it is the reciprocal of this function which must be integrated in the Picard iteration procedure. A means of avoiding this problem is to use a simple function such as a polynomial or an exponential function to provide simple upper and lower bounds to the function 1/El(x) in a given range. If these bounds are close, the bounds constructed on the simpler problem by the Picard iteration procedure will converge to close bounds on the actual problem. In the above analysis the upper and lower bounds on the transfer matrix elements for a static or dynamic beam are constructed independently. When elastic or rigid supports are included as a part of the structure, these bounds become coupled, Coupling of upper and lower bounds can occur for the beam alone if negative elements occur in the [A] matrix : e.g., for a beam-column analysis. From the transfer matrix bounds, one can proceed to the development of both upper and lower bounds to the entire state vector. Bounds on all
BOUNDSFOR VARYINGGEOMETRYBEAMS
447
quantities of the state vector can be obtained to essentially the same degree of precision. It is shown that the actual time of computation for generating the upper and lower bounds is significantly less than that required for an ordinary transfer matrix analysis using the RungeKutta procedure. The accuracy of the solution is apparent when upper and lower bounds are available, whereas this is not the case for the numerical solutions. It is necessary in the ordinary transfer matrix approach to use techniques such as the delta matrix method [2], formulation of a super matrix equation [9], etc., to circumvent numerical difficulties associated with the analysis of large and multispan structures. Such techniques can also be employed when bounds on the transfer matrix elements are available. It is important to use all possible information that the differential inequalities can supply in order that the above techniques will give bounds which are sufficiently accurate to be of use to engineers. ACKNOWLEDGMENTS This research was sponsored by the National Science Foundation under Grant No. GK-40589 and by the Iowa State University Engineering Research Institute. Thanks are due to Mr P. V. Thangam
Babu and Miss Cecelia Byers for assistance
with this manuscript.
REFERENCES 1. T. J. MCDANIEL and J. P. HENDERSON1974 The Shock and Vibrution Digest 6-l (Part I), -2 (Part II). Review of transfer matrix vibration analysis of skin stringer structure. 2. E. C. PESTELand F. A. LECKIE1963 Matrix Methods in Elastomechanics. New York: McGrawHill Book Company, Inc. 3. Y. K. LIN 1967 Probabilistic Theory of Structural Dynamics. New York: McGraw-Hill Book Company, Inc. 4. Y. K. LIN and T. J. MCDANIEL1969 Journal of Engineering for Industry 91,1133-l 141. Dynamics of beam-type periodic structures. 5. C. A. MERCERand C. SEAVY1967 Journal of Sound and Vibration 6,149-l 62. Prediction of natural frequencies and normal modes of skin-stringer panel rows. 6. D. J. MEAD 1971 Journal of Engineering for Industry 93, 783-792. Vibration response and wave propagation in periodic structures. 7. D. J. MEAD and G. SENGUPTA 1970 USAF Report AFML-TR-70-13. Propagation of flexural waves in infinite, damped rib-skin structures. 8. J. P. HENDERSON1972 Ph.D. Thesis, Ohio State University. Vibration analysis of curved skinstringer structures having tuned elastomeric dampers. 9. J. P. HENDERSONand T. J. MCDANIEL 1971 Journal of Sound and Vibration 18, 203-219. The analysis of curved multi-span structures. 10. T. J. MCDANIEL 1971 Journal of Aircraft 8, 143-149. Dynamics of circular periodic structures. 11. T. J. MCDANIEL 1792 Journal of Sound and Vibration 23, 217-227. Dynamics of non-circular stiffened cylindrical shells. 12. W. F. HUNTER 1970 NASA TN-D 6064. Integrating matrix method for determining the natural vibration characteristics of propeller blades. 13 T. J. MCDANIEL and J. D. LOGAN 1971 Journal of Sound and Vibration 19, 39-48. Dynamics of cylindrical shells with variable curvature. 14. W. WALTER1970 Diflerential and Integral Inequalities. New York: Springer Verlag. 15. M. H. PaorrEa and H. S. WEINBERGER1967 Maximum Principles in D@rential Equcttions. New Jersey: Prentice-Hall, Inc. 16. V. LAKSHMIKANTHAM and S. LEELA1969 D@erentialandIntegral Inequalities, Volume I, Ordinary Differential Equations. New York and London: Academic Press. 17. B. FRAEJIS DE VEUBEKE 1965 Proceedings of Conference on Matrix Methods in Structural Mechanics, AFFDL-TR-66-80, 863-886. Bending, and stretching of plates-special models for upper and lower bounds. 18. N. W. BAZLEY, D. W. Fox and J. T. STADTER1967 Zeitschrift fur angewandte Mathematik und Mechanik 47, 191-198. Upper and lower bounds for the frequencies of rectangular clamped plates.
448
T. I. MCDANIELAND V. R. MURTHY
19. K. NICKEL 1962 Ingenieur Archiv 31, 85-l 00. Eine einfache Abschatzung fur Grenzschichten. 20. E. ADAMS 1967 Acta Mechanica 3, 268-277. Rigorous bounds of exact solutions of boundary layer problems of Prandtl type. 21. F. JOHN 1965 Courant Institute of Mathematical Sciences, New York University, Research Report No. IMM-NYU 308. Estimates of the error in the equations of non-linear plate theory. 22. J. R. CANNON 1966 Proceedings of an Advanced Symposium Conducted by the Mathematics Research Centre, United States Army, at the University of Wisconsin, Madison. 21-54. Numerical solutions of nonlinear differential equations. 23. C. A. BELL and F. C. APPL 1974 JournalofAppliedMechanics 40, 1097-l 102. Bounds for initial value problems. 24. E. T. CRANCH and A. A. ALDER 1955 Journal of Applied Mechanics 21, 103-108. Bending vibrations of variab!e section beams.