- Email: [email protected]
Thermally induced vibrations of beam structures
COMPUTER METHODS IN APPLIED @ NORTH-HOLLAND PUBLISHING
MECHANICS COMPANY
AND ENGINEERING
21 (1980) 337-35.5
THERMALLY INDUCED VIBRATIONS OF BEAM STRUCTURES G.D. MANOLIS* Department
and D.E. BESKOS**
of Civil and Mineral Engineering, University Minnesota 55455, USA Received Revised
19 January 9 March
of Minnesota,
Minneapolis,
1979 1979
A general numerical method is developed for determining the dynamic response of beam structures to rapidly applied thermal loads. The method consists of formulating and solving the dynamic problem iq the Laplace transform domain with the aid of dynamic stiffness influence coefficients defined for a beam element in that domain and of obtaining the response by a numerical inversion of the transformed solution. Thus, the solution of the associated heat conduction problem, usually obtained by Laplace transform and needed for the computation of the thermal load, can be used in its transformed form. The effects of damping and of axial compressive forces on the structural response are also studied. Three examples are presented in detail to illustrate the proposed method and demonstrate its advantages.
1. Introduction The conventional thermal stress analysis of structures ignores the effect of inertia and, by using a predetermined time-dependent temperature distribution, essentially deals with a static problem where time is reduced to the role of a parameter. For slender structures under rapid surface heating however, such as beams, plates and shells, the inertia effects become significant and lead to thermally induced structural vibrations. Boley [l], [2] was the first to take structural inertia into account and analytically treat the problems of thermally induced vibrations of simply supported elastic beams and plates. His work was followed by the analytical studies of Kraus [3] on thermally induced vibrations of thin elastic spherical shells and of Yu [4] and Frisch [5] on observed thermally induced vibrations of flexible satellite booms subjected to variable heating. Boley [6], reconsidering the problem, derived a simple approximate formula for the ratio of the maximum dynamic to static deflection of heated beams and plates, and he then generalized it to include the effects of damping and of axial or in-plane loads. In addition to the aforementioned thermally induced vibrations in aerospace structures due to surface heating, similar vibrations due to radiation heating have been observed and studied in nuclear structures [7], [8]. A review article by Bargmann [9] discusses both types of thermally induced vibrations. This paper is concerned with thermally induced vibrations of structures composed of beams which are exposed to rapid surface heating. For complex beam structures, such as continuous beams and gridworks or frameworks, the analytical approach [l]-[5] based on the MindlinGoodman method [lo] for solving equations with time-dependent boundary conditions is not *Graduate Student. **Associate
Professor.
338
G.D. Manolis, D.E. Beskos, ~e~u~ly
induced oibrutio~ of beam structures
efficient, and numerical procedures should be used. In this work a general numerical method of analysis is developed which is capable of determining the elastic dynamic response of any plane beam structure to any time-dependent thermal loading. The method consists of applying the Laplace transform with respect to time to the governing equations of bending motion of an elastic beam element, solving for the transformed displacement function and constructing stiffness coefficients in the transformed domain. The problem is then formulated and solved in the transformed domain, and the transformed solution is inverted numerically to yield the dynamic response. The transformed thermal loading is usually obtained by solving the corresponding heat conduction problem in the Laplace transform domain. This method of solution of dynamic problems was first presented by Beskos and Boley [ll] and more recently generalized and extensively treated by Narayanan and Beskos [12]. The most important problem in this method is the numerical inversion of the Laplace transform, which is discussed in detail [12]. Here the method based on the trigonometric expansion as described by Papoulis [13] and the method based on sine and cosine transforms as described by Durbin El43 are utilized, The second method is more time-consuming than the first one but is more accurate for short, as well as long, time solutions. Other effects such as damping, rotatory inertia and shear deformation, as well as the effect of axial loading on flexure, can be easily considered by constr.ucting appropriate stiffness coefficients in the Laplace transform domain to include them. Here the effects of damping and of axial loads on the response are studied. Illustrative structural examples of beams and gridworks are presented to demonstrate the use and importance of the method, which presents distinct advantages over other well-known methods as far as accuracy, simplicity and efficiency are concerned.
2. General formulation and solution This section describes in detail the proposed numerical method for formulating and solving the general problem of thermally induced vibrations of beam structures exposed to rapid surface heating. Consider a uniform linear elastic Bernoulli-Euler beam element l-2 of length L, bending rigidity EI, cross-sectional area A and mass per unit of length m under rapid surface heating as shown in fig. 1. On the assumption that the resulting thermal loading is concentrated at nodes 1 and 2, the equation of flexural motion for the beam element l-2 is given by
where z, = u(x, t) stands for the lateral beam deflection along the positive y direction of fig. 1, and t denotes time.
Fig. 1. The flexural
beam element
in the Laplace
transform
domain.
G.L>. ~anol~s, D.E. Beskos, penalty
By following [Ill the Laplace transform initial conditions to yield
undyed vibrationsof beam secures
with respect to time is applied to (1) under zero
$+4zc4u‘=o, where a(x, S) =
I0
339
(2)
m U(X,t) e-“’ di
(3)
and 4K4 = ms2/EI. The
(4)
solution of eq. (2) is given in [15] as 6(x, S) = eKX[A COSKX) + B sin(&)] + eeKX[C cos(Kx) + D sin@)],
(5)
where K is given in terms of s by (4), and A, B, C and D are constants of integration. The beam element 1-2 has two degrees of freedom per node with positive directions for them as well as for their co~esponding nodal forces as shown in fig. 1. Thus, using 6(x, s) as the displacement function in a stiffness formulation, the beam element stiffness matrix can be constructed in the Laplace transform domain following standard procedures. The following relations are utilized in constructing the element stiffness matrix: 8(x, S) = -d;iildx, V(X, S) = -EI(d3E/dx3),
(6)
ti(x, s) = EI(d2$dx2), where 8 is the rotation, and v and $i are the shear force and the bending moment, respectively, with positive directions taken as in f16]. Finally the dynamic stiffness equation for the beam element l-2 in the Laplace transform domain takes the following static-like form:
(7)
In eq. (7), where the subscript T indicates thermal load, the dynamic stiffness influence coefficients in the Laplace transform domain are given by the following expressions: Dl, = 2NK2(y2 + 4ysc - l), 62, = Mq-y* d,,
+ 2y(l-
2s2) - 11,
= (4IvK”/vg[-y’(s3
lY,l = (2NK/%5)(s3
+ SC2+ c) - y(s3 + SC2- c)],
+ SC2+ s)(y - y2),
L522= N(y2 - 4ysc - l), n+* = (2rJ;tv7)[y’(s’
+ SC2- c) + y(s3 + SC2+ c)],
Is,, = -&,
D33 = Ds,,,
D43 = -LTZ,,
6,
= D2*,
(8)
340
C.D. Man~~is, D.E. Beskos, ~er~a~~y
induced vibrations of beam structures
where N = 2J51K/[y2 - 2y(1+ 2s2)+ 11, c = cos(KL),
s = sin(KL),
y = exp(2KL).
(9
With the introduction of the coefficients rS, the dynamic problem of a plane beam structure in flexural motion due to thermal loading is expressed in the Laplace domain in the static-like form F(S)T = ~(~)~(~),
(10)
where F(s)T and G(S) are the Laplace transformed thermal load and displacement vectors, respectively, while B(s) is the transformed dynamic stiffness matrix of the structure obtained by an appropriate superposition of element stiffness matrices as the one in (7). After application of the transformed boundary conditions, the transformed solution G(S) of eq. (10) is obtained, for a sequence of values of S, by a standard matrix inversion procedure usually performed numerically. From the above solution the response u(t) is obtained, for a sequence of values of t, by a numerical inversion of the Laplace transformed displacements. Heating of a beam element produces a thermal nodal loading consisting in general of thermal bending moments and axial forces. Thermal axial forces can be disregarded on the assumption that axial deformations are not restrained, while thermal moments are given by MT=Ea
Ty W
(11)
A
where (Yis the coefficient of thermal expansion, which is assumed to be constant, and T is the temperature dist~bution in the beam, in general a function of X, y and t. The temperature dist~bution T is obtained by solving the heat conduction problem defined for the heated structure. Consider an isotropic solid body with a constant thermal diffusivity U, mass density p and specific heat c, subjected to arbitrary thermal boundary and initial conditions on its surface and with internal heat generation 4 per unit time per unit volume, which is of importance in problems of radiation heating [7], [8] but is taken as zero for problems of surface heating as in the present case. The general equation of heat conduction is (see [17])
Uv*T+-g, PC
(12)
where T is the temperature of the body, varying with space and time, and V2 is the Laplace’s operator_ Eq. (12) is supplemented with appropriate boundary and initial conditions to form a well-posed boundary value problem. The five principal boundary conditions of heat conduction are described in detail in [17]. For the problem of surface heating of a uniform flexural beam element of rectangular cross-section with depth h and width b one has that 4 = 0 and T = T(y, t) in (12) so that this equation becomes (13)
G.D. ~a~olis,
D.E. Beskos, verbally
induced uibraiio~s of beam sfruc~ures
341
If 7’ is a function of x as well, one can afways subdivide the structure into elements with constant temperature along the beam element length. In the present problem the initial condition is assumed to be
0) = 0,
VY,
(14)
while the boundary conditions K
are
5 =-Q(t)
at
KO
y = h/2,
at y =
ay-
(15)
-h/2,
indicating application of a heat input Q(t) on the top beam fibers and insulation of the bottom ones. The constant K in (15) is the thermal conductivity of the material and is related to u by
One usual way of solving the heat conduction boundary value problem consisting of (13)-(H) is by the method of Laplace transform. With the application of the Laplace transform with respect to time eqs. (13~(15) are transformed into the system (d2ndy2) - (s/u)l?; = 0, K(dmdy) = -b(S),
at
y = h/2,
dp,‘dy = 0
at
y = -h/2,
F = F(y,
s) = fin Tfy,
t) e-“’ dt,
&)
= lrn Q(t) e-“’ dt,
(17)
0
which yields the following solution for the temperature F(y, s) = A eAy+ B eeAy,
in the Laplace transform domain: (18)
where h = %5/U, A = -(~/~~)~e~~~~/(lB = Since the transform by Laplace of bending
-(@hK)[ehh’2/(1
-
eU”)],
(19)
e2Ah)].
problem of the structural vibration is formulated and solved in the Laplace domain, there is no need to obtain the temperature distribution in the time domain inversion. In view of (ll), the nodal thermal loading in the Laplace domain consists moments which for a rectangular beam are of the form M;_(s) = =b
1-y; ?‘(Y, 4~ dy,
where F(y, s) is given by (18) and (19).
(20)
342
G.D. Manolis,
D.E. Beskos,
Thermally
induced vibrations of beam structures
The heat input Q(t) can be any fast-rising function of time such as or a gradually applied (ramp-type) one with small rising time. In nonlinear function of time which can be approximated by a piecewise case, of course, it is preferable to compute o(s) and consequently following [ 121.
a suddenly applied one fact, Q(t) can be any linear function. In that MT(s) numerically by
3. Numerical inversion of Laplace transform Let f(t) be a real function of t, with f(t) = 0 for t -=L 0. The Laplace transform f(s) of f(t) is defined by
T(s)= 6 f(t)
e-“’ dt,
(21)
and its inversion formula is given by
where p > 0 is arbitrary but is greater than the real part of any of the singularities of f(s), and s is a complex number with Re(s) I j3 > 0. In many engineering applications f(s) is either too complicated to be inverted analytically by evaluating the integral of (22) or is given in numerical form. In those cases a numerical inversion of the Laplace transform is imperative. Reference [12] discusses this problem in detail and investigates the accuracy, computational efficiency and range of applicability of eight methods of numerical Laplace inversion as applied to dynamic problems. Among those methods two are briefly presented and used in this work. The first method is due to Papoulis [13] and determines f(t) in terms of the values of f(s) on a finite sequence of equidistant points s, on the real axis by the formula f(o) = 5 C, sin(2n + l)O, n=O
(23)
where N is a finite positive integer number, the coefficients C,, are given in [13] as functions of the values of T(s)), and 0 = cos-‘(e-“‘),
(24)
where S > 0 is a constant selected by the user. This method has the advantage that it works with real data and provides easy-to-use expressions for the evaluation of the C”, thus requiring small amounts of computer time. Unfortunately, the accuracy of the method is only acceptable for early times, and improvement of its accuracy by increasing hi is not possible because the computation of the C, becomes an unstable process for large N [18]. The second method of inversion is due to Durbin [14] and is actually an efficient improvement
G.D. Manofis, D.E. Beskos, Thermally induced vibrations of beam structures
343
of the method of Dubner and Abate [19], which is based on the finite Fourier cosine transform. Durbin [14] combined both finite Fourier cosine and sine transforms to obtain the inversion formula f( tj )
x
eaiAr [-+Rev@)]+Re{y
(A(n)+iR(n))W”]], tl=O
where
Re{f@
+ i(n + 1N) F
B(n) = 5 Im{f(P I=0
+t i(n + fN)+
A(n) = i
1, tz = 0, 1, . . . )
and f(t) is computed for N equidistant points jAt = jT/N, j = 0, 1,. . . , N - 1. It is suggested that for L x N ranging from 50 to 5000 one should select fit = 5 to 10 for good results, where t^ is the total time interval of interest. The computations involved in (25) are performed by employing the fast Fourier transform algorithm of Cooley and Tukey [20]. The accuracy of Durbin’s method [14] is very high even for late times, and its only disadvantage is the fact that it works with complex data and thus is more time-consuming than the method of Papoulis [13]. This disadvantage is however offset by its high accuracy, as the numerical examples in this work will clearly demonstrate.
4. Damping and axial force effects The effect of damping on the thermally induced motion of a beam has been investigated by Yu f4] and Boley [6]. Yu [4] considered internal viscoelastic damping of the beam material by employing the constitutive law of the standard linear solid and installing a viscous-fluid damper at the tip of the cantilever beam when modeling a satellite boom. Boley [6] studied the damping effect by considering viscous damping in a single degree of freedom system approximating beams or plates. In this work both internal viscoelastic damping and external viscous damping are considered. Internal viscoelastic damping is accounted for by assuming that the beam material obeys the constitutive law of a Kelvin viscoelastic model: CT= E(E + j-k),
(27)
where u is the stress, E is the strain, g is the strain rate, E is the moduius of elasticity, and f is the internal damping coefficient of the material (which is usually very small for metals). Application of the Laplace transform as defined by (21) under zero initial conditions on (27) yields (T = E(l + fS)Z,
(28)
G.D. Manolis, D.E. Beskos. Thermally induced vibrations of beam structures
344
where Cr and E are the transformed stress and strain. Eq. (28) indicates that internal viscoelastic damping can be very easily taken into account if E is replaced by E(l + fs) in eqs. (9) and (20). In the case where there is external viscous damping eq. (1) is replaced by (29) where c is the viscous damping coefficient. Application defined by (3) yields
of the Laplace transform
to (29) as
$$+4K4C=0,
(30)
4K4 = (ms’+ cs)/EI.
(31)
where
Thus, in view of the fact that (2) and (30) are identical in form, external viscous damping can also be very easily taken into account if (4) is replaced by (31). The effect of axial loads or beam-column effects on the thermally induced motion of a beam has been approximately treated by Boley [6], who utilized the fact that, for simply supported beams under certain types of loads, the magnification of the static or dynamic deflection under an axial compressive load P can be approximated by the formula 1
& =1 - (PIP,,)’
(32)
where PC, is the buckling load. In this work the effect of axial loads on the response of beam structures to thermal loading is determined in an exact way by constructing the coefficients Dii which include the effect of axial loads on flexure. In the presence of a constant compressive force P, eq. (1) is replaced by EI ~+P$+m$=O, a4v which, after application
(33)
of the Laplace transform as defined by (3) becomes (34)
where 2~ = PIEI,
v = ms2/EI.
(35)
The solution of (34) is given by a(~,
s) = A e*‘”+ B e’2x+ C eAp + D e”“,
(36)
G.D. Manolis, D.E. Beskos, Thermally induced vibrationsof beam structures
where A, I?, C, D are constants numbers given as
of integration,
345
and AI, AZ, hJ, A4 are in general complex
A,=&-,
A2 = -J-CL
+ e,
A3=+p-m,
A4 = -d-p
-m.
(37)
On the basis of the displacement function (36) the D, including the effect of P can be constructed as in section 2. Here, however, due to the lengthy algebraic computations involved, explicit expressions for the Dij are not derived. Only the constants of integration A, B, C and D of (36) are explicitly determined and given in the appendix for the necessary displacement states which define the Dij. For a beam-column element the Dij are computed numerically on the computer with the aid of (36) and v(x, s) = -EI(d3G/dx3) - P(dG/dx), n;i(x, s) = EI(d2G/dx2).
(38)
The case of a tensile axial force P can be very easily taken into account by simply replacing P by -P in all the previous expressions. The combined damping and beam-column effects can also be treated on the basis of eqs. (33)-(38) either by replacing E by E(l + fs) for the case of internal viscoelastic damping or by defining v as v =
(ms2+ cs)/EI
(39)
for the case of external viscous damping.
5. Structural examples This section describes the detailed solution of three structural examples in order to illustrate the use of the proposed numerical method, demonstrate its advantages and provide qualitative as well as quantitative results for the problem of thermally induced vibrations of plane beam structures. Example 1 Consider a uniform simply supported beam of length L, bending rigidity EI, mass per unit of length m and of rectangular cross-section with width b and depth h as shown in fig. 2. A
Fig. 2. A simply supported
beam under suddenly
applied
heat input.
346
G.D. Manolis,
D.E. Beskos.
Thermally
induced
vibrations of beam structures
step heat input Q(t) of intensity QO, uniformly distributed along the length of the beam, is assumed to be applied along the edge y = h/2, while the edge y = -h/2 is insulated, as fig. 2 shows. The deflection and bending moment of this beam have been analytically obtained by Boley [l], who presented his results in terms of the following dimensionless parameters:
T =
utlh2,
VW
V(& T) = 7r4~v/(192Q~L2), mT(& T) =
(4Oc)
n4~hfT/(192EI@),
(404
where z, and MT stand for deflection and thermal moment, respectively, and all the other quantities which appear in (40) have been defined in section 2. It was shown in [l] that the governing nondimensionless parameter of the problem, which is an index of the importance of the inertia effect, is B. This is defined by (40a) as the square root of the ratio of the characteristic thermal time h2/u to the characteristic mechanical time (-the first natural period of vibration). It was found in [l] that the importance of the inertia effects increases with decreasing B so that for B = 0 they actually prevent any beam deflection. Also, as B approaches infinity, the inertia forces disappear and the static solution results. Furthermore, it was also shown in [l] by a quantitative analysis on an aluminum beam that small values of B for which inertia effects are important correspond to very thin beams. In this work, on the assumption that the beam material is aluminum, the following numerical data was adopted for the computation of the midspan deflection and bending moment of the beam l-3 of fig. 2: b =
1 in,
h =
0.1539 in
L = 40 in,
m = 4 x 10d5 lb sec2/in2,
E = 10.66 X lo6 lb/in2,
u = 0.1333 in21sec,
K
cy = 12.25 X 1O-6l/“F,
Q. = 1 Btu/in2 sec.
=
0.0027 Btu/in set “F, (41)
On the basis of the above data one has that B= 1,
T = 5.62789&
V = O.O699v,
mT = O.O34966M,
(42)
In order to determine the midspan dynamic deflection n2 of the beam l-3 of fig. 2, a node is introduced at midspan, thus dividing the beam into two elements. Following the general procedure of section 2 and utilizing the boundary conditions & = iT3= 0, 0, = 0 and 0, = -&, the deflection C2 due to nodal thermal loads M1, = -Z&, = i& and l&- = 0 takes the form
G.D. Manolis, D.E. Beskos, Thermally induced vibrations of beam structures
347
where the Oij are given by (8), and A& is given by (18)-(20) for a(s) = Qols. The response ~)~(t)was obtained by numerically inverting the Laplace transformed solution U2given by (43). The methods of Papoulis [13] and Durbin [14], as described in section 3, were both employed, and the results are shown in fig. 3 together with the exact analytical results of [l]. It is apparent in fig. 3 that, while Durbin’s results indicate a perfect agreement with the exact solution, this is not the case with Papoulis’ results, which present a reasonable agreement only for the first cycle of vibration. The CP computer time spent on the CDC Cyber 74 computer of the University of Minnesota Computer Center was 2.118 set for Durbin’s and 1.285 set for Papoulis’ methods, respectively, for 250 time points. The values of 6 = 5.866 and N = 9 for Papoulis’ method and L = 1, N = 250, t^= 1.0 set and pt^ = 6 for Durbin’s method were used in this example. Due to the high accuracy results of Durbin’s method all further numerical Laplace inversions were performed by that method, even though more computer time is required by that method _than by Papoulis’ method. The bending moment at the midspan of the beam l-3 of fig. 2 was also computed by observing that
M,=M;i-AT;,
(44)
where the dynamic moment A?,” is obtained by considering the beam element l-2 and solving for the internal element forces in terms of known nodal displacements and has the form
and the static moment A?; is equal to fi, where A& is given by (18)-(20) for a(s) = QJs. The bending moment A& given by (44) was inverted numerically by Durbin’s method, and the results show excellent agreement with the exact ones of [l], as fig. 4 demonstrates in a dimensionless plot of m2 = 7r4~M2/192QcuEI versus 7 = ut/h*.
-
Bolcy PI
--*-*
Papoulls
0000
Durbln
[Ia] [I41
Fig. 3. Dynamic midspan deflection of the simply supported beam for B = 1, f
= 0.
348
G.D. Manolis, D.E. Beskos, Thermally induced vibrations of beam structures
m
-
a3
B&y
* * 0 The
[I] method
02
Fig. 4. Dynamic midspan bending moment of the simply supported beam for I? = 1, ,f = 0.
The effect of internal viscoelastic damping on the midspan deflection is portrayed by fig. 5. The computations were done on the basis of (43), with E being replaced by E(l +fs) everywhere. Because the internal damping coefficient f is very small for metals and ranges between 0.005 and 0.020 (see [21]), the values 0.005, 0.010 and 0.015 were used in this example. Fig. 5 shows clearly the attenuating effect of damping on the dynamic response, which increases rather rapidly with increasing values of f and time T. Eventually, after some time the damped dynamic deflection approaches the static one which corresponds to the case of zero inertia effects. This indicates that damping is in general an effective way to reduce inertia effects in a thermal stress structural problem. 4 V
f
0.06 -
0
0.5
1.0
1.5
-c
Fig. 5. Dynamic midspan deflection of the simply supported beam for various values of the intertia index B and the damping f.
G.D. Manolis, D.E. Beskos, Thermally induced vibrations of beam structures
349
A
v,
(If-l)
3.0 -
This
method
2.0 -
1.0 -
c
0
Fig. 6. Dynamic
midspan
deflection
0.05
0.20 t
0.15
0.10
of the simply supported
(Set)
beam for B = 1, f = 0 and various
values of P/P,,.
Finally fig. 6 depicts the amplification of the midspan deflection v2 of the beam of this example due to the presence of an axial compressive force P
+
t
Q(t) QO
Q(t)=0
\/
1
i/l
1 4 l_ A
El m T
+-C--+-C_, 2 Fig. 7. A two-span
continuous
El m y 2
beam under suddenly
applied
heat input.
G.D. Manolis, D.E. Beskos, Thermally induced vibrations of beam structures
350
obtain the following expressions for the rotations at the supports of the continuous
beam:
-0, = -Ii&& &=o,
(46)
where fir is given by (18j(20) for o(s) = Qofs, and & is given by (8). A numerical Laplace inversion of 6, yields o,(t), which is plotted in fig. 8 for f = 0 and f = 0.015. Fig. 8 shows that the inertia effect is not so significant in this case as B,(t) oscillates with small amplitude about its static value. In fact, for a damping value of f = 0.015 there is virtually no oscillation, and o,(t) rather rapidly approaches the static solution.
Fig. 8. Dynamic rotation @x(t)of the continuous beam for B = 1, f = 0 and 1.5%.
Example 3 Consider a plane horizontal square gridwork simply supported at its four corners as shown in fig. 9. The gridwork consists of beam elements of length L, = 24 in and depth h = 0.2 in and with all other geometrical and material properties as in (41). It is assumed that a suddenly applied heat input acts on the top fiber of all the beams of the gridwork, while their bottom fibers are insulated. It is further assumed that axial and torsional deformation of the beams can be neglected so that element equilibrium equations in the Laplace transform domain are of the form (7). In this example, due to the large number (23) of free degrees of freedom of the gridwork, the formulation and solution of the problem was done exclusively by computer. A general computer program was written in FORTRAN capable of performing a complete dynamic analysis of a gridwork under thermal loading on the basis of the proposed method as described in section 2. The program constructs the total dynamic structural stiffness matrix fi and solves equation (10) numerically by Gauss-Jordan elimination for a sequence of complex values of s. The
G.D. Manolis, D.E. Beskos, Thermally induced vibrationsof beam structures
351
Fig. 9. A plane horizontal square gridwork under suddenly applied heat input.
transformed solution thus obtained in discrete form is inverted numerically by Durbin’s method to yield the dynamic response. Using this program, the vertical deflection at node 9 of the gridwork of fig. 9 was determined as a function of time, and the result is shown in fig. 10. The total CP computer time for this computation was 20 sec.
20-
1.0-
i
00
0.1
0.2
0.3
t mx)
Fig. 10. Dynamic vertical central deflection of the gridwork for B = 1, f = 0.
6. Concl~ions On the basis of the preceding discussion one can draw the following conclusions: 6.1
A general numerical method for determining the thermally induced dynamic response of elastic beam structures has been developed. The method formulates the dynamic problem in the Laplace transform domain-in a static-like form by using Laplace transformed dynamic stiffness influence coefficients -ct, and solves it (there) numerically. The dynamic response is finally obtained by a numerical inversion of the transformed solution.
352
G.D. Manolis, D.E. Beskos, Thermally induced vibrations of beam structures
6..2 Only Laplace transformed dynamic thermal loads are used in this method. Thus the associated transient heat conduction problem under any kind of boundary conditions has to be formulated and solved in the Laplace domain only, and there is no need for its time domain solution. 6.3 In the framework of the linear theories employed here the use of the Dij coefficients leads to the “exact” solution of the dynamic problem since the displacement function used for their construction is the exact solution of the transformed equation of motion. Thus the present method can provide a basis for comparing the accuracy of other approximate methods such as the finite element method or the method of Boley [6], which is based on the assumption that the first-mode response is legitimate for many practical problems. 6.4 Various effects on the response of beam structures such as damping, axial forces, axial or torsional deformations, rotatory inertia and shear deformation can be very easily taken into account by constructing appropriate Dij coefficients which include these effects. In this work the effects of damping and axial forces were considered and studied. As far as damping is concerned, use of Dij coefficients permits consideration of different amounts of damping in different beam elements, thus enabling one to effectively control the response by appropriate damping changes in particular structural elements. 6.5 The proposed method is quite general and capable of solving any dynamic problem of beam structures and not just thermally induced vibrations. From this point of view it can be considered as a generalization of the method proposed by Beskos and Boley [ll]. Because in problems of thermally induced vibrations the principal results for beams and plates are almost identical regardless of the aspect ratio of the plate or the rapidity of the heat application and even of the support conditions [6], the results obtained by the proposed method for beam structures can be quite informative about the response of their associated plate structures. Apart from that, the basic ideas of the method can be extended to treat dynamic problems of plate and shell structures directly. 6.6 The proposed method appears to be better than the conventional finite element method in conjunction with either modal analysis or numerical integration because i) it provides the exact solution of the dynamic problem and not an approximate one as in the finite element method (which discretizes the continuously distributed mass of physical beam elements), and ii) it considerably simplifies the solution of the associated heat conduction problem. Furthermore, the method does not require prior knowledge of natural frequencies and modal shapes as in modal analysis. The only disadvantage of the method has to do with the numerical Laplace
G.D. Manolis, D.E. Beskos, Thermally induced vibrations
of beam
structures
353
inversion, where one has to decide between relatively low accuracy with less computer time and high accuracy with more computer time and select his inversion algorithm on this basis. Reference [12] indicates that, for a given dynamic input, numerical integration techniques such as Newmark’s p-method or Wilson’s O-method require in general less computer time for a dynamic problem than numerical Laplace transform techniques for the same degree of accuracy. The latter techniques however have the characteristic properties that they permit one to perform numerical inversion for just one nodal displacement of interest independently of the others, or determine a particular nodal displacement at a desired time without prior knowledge of its values at all previous times. From these view points numerical Laplace transform techniques are more efficient than numerical integration methods. Also, as [22] points out, Laplace transform techniques based on algorithms of numerical inversion working with real data seem to be more economical than numerical integration for large systems to be studied over a large time domain because no more than 20 system resolutions need to be computed.
Acknowledgement
The authors are grateful to the University of Minnesota Computer Center for making its facilities available to them and to Dr. F. Durbin for supplying them with his numerical Laplace inversion computer program.
Appendix. Determination displacement states
of the constants of integration A, B, C and D of (36) for the various
1st displacement state. 6(0, s) = 1, ~‘(0, S) = I?(& s) = IT’& S) = 0, where the prime indicates differentiation with respect to x. For this case the use of (36) results in the system A+B+C+D=l, hlA + A,B + h,C + AbD = 0, ehlLA + eAzLB + eA*C + eALD = 0, Al
eAILA + A2 ehzLB + A3 eA*C + A, e”&D = 0,
(A4
which has the solution A = (l/A)[hz(A4 - A3)e(*3+*4)L + A3(A2- A4)e(h2+A4)L + A4(A3- AZ)e(“2+A3’L], B = (l/A)[A1(A3 - A4)e(h3+A4)L + A3(A4- A,) e(h1CA4K + A4(Al- AS)eQ’+*3)L], C = (l/A)[AI(A, - AZ)e(A2+h4)L + A2(Al- A4)e(h1+*4)L + A4(Az- AI) e(A*+*2)L], D = (l/A)[A,(A,-
As) e(A2+*3w + A,(A, - A,) e(A1+*3)L + Ag(AI- AZ)e(h’+*2)L],
64.2)
where A = sAlA + (A, - A3)(A4- AZ)[e(AlfA3ZL + e(A2+r4)L] + (A, - A4)(A2- A3)[e(A1+h42L + e(A2+*3)L].
64.3)
354
ti.D. Manolis, D.E. Beskos, Thermally induced vibrations of beam structures
2nd displacement state. a(O, s) = 0, U’(0, s) = 1, U(L, S) = c’(L, S) = 0. For this case the use of (36) results in the system A+B+C+D=O, A,A + h2B + h,C+ h,D = 1, eAILA + eAzLB+ eAXLC+ e*&D = 0, A,
e”lLA + A2 eA2LB + A3eA’LC +
Aq e**D
(A.4
= 0,
which has the solution
A =
(l/A)[(A,
-
B
=
(l/A)[(A, - A3)e(*3+*4jL + (A, - A4)e(*l+*d)L + (A3_ A,) e(*l+~T)~],
C
=
(l/d)[(hz-
D
=
(l/A)[(A, - A*)e(A*+AX)L + (A, - AX)e(*l+*J)L+ (A, _ A,) e(*~+~~)~],
Ad)
e(A3+A4)L + (Ad-- hZ)e(*2+*4)L
+ (AZ-
~~)~f*z+AdL],
Aq)e(A2+A4)L + (A4- A,) e(A1+*4)L + (A, _ AZ)e(*~+~~)~~,
(A.3
where A is again given by (A.3). Eqs. (38) computed at nodes 1 and 2 of the beam element of fig. 1 become V(O,s)=-EI(AA:+BA;+CA:+DA:)-P(AA1+BA,+CA,+DA,), M(o,s)=EI(AA:+BA;+~A:+DA:), v(L, s) = -EI(AA
f e*IL +BA~eA2L+CA~eA~+DA~eA4L)
- P(AA, eAIL+ BAZeAzL+ CA, e*& + DA, e*&), I@(& s) = EI(AA: eAIL+ BA$ eA2L+ CA: e*$ + DA: e*&).
(A.6)
Using (A.6), the Oij coefficients for the 1st and 2nd displacement states can be obtained by utilizing the values of the constants A, B, C and D for the corresponding displacement state as given by (A.2) and (AS). Thus the 1st ,and 2nd displacement states in conjunction with the proper sign convention consideration provide the first and second columns of the beam element stiffness matrix 0. The third and fourth columns can be very easily obtained from the first two by using the Betti-Maxwell reciprocal relations. References [l] [2] [3] [4]
B.A. Boley, Thermally induced vibrations of beams, J. Aerospace Scs. 23 (1956) 179-181. B.A. Boley and A.D. Barber, Dynamic response of beams and plates, J. Appl. Mech. 24 (1957) 41Wl16. H. Kraus, Thermally induced vibrations of thin nonshallow spherical shells, AIAA J. 4 (1966) 500-505. Y.Y. Yu, Thermally induced vibration and flutter of a flexible boom, J. Spacecraft and Rockets 6 (1969) 902-910. [S] H.P. Frisch, Thermally induced vibrations of long thin-walled cylinders of open section, J. Spacecraft and Rockets 7 (1970) 897-905.
G.D. Manolis, D.E. Beskos, Thermally induced vibrations of beam structures [6] B.A. Boley, [7] [8] [9] [lo] [ll] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]
355
Approximate analyses of thermally induced vibrations of beams and plates, J. Appl. Mech. 39 (1972) 212-216. W.C. Lyons, Comments on heat induced vibrations of elastic beams, plates, and shells, AIAA J. 4 (1966) 1502-1503. H. Bargmann, Dynamic thermal shock resistance, pp. 174-181, in: J.L. Zeman and F. Ziegler (Eds.), Topics in applied continuum mechanics (Springer, Berlin, 1974). H. Bargmann, Recent developments in the field of thermally induced waves and vibrations, Nut. Eng. Design 27 (1974) 372-385. R.D. Mindlin and L.E. Goodman, Beam vibrations with time-dependent boundary conditions, J. Appl. Mech. 17 (1950) 377-380. D.E. Beskos and B.A. Boley, Use of dynamic influence coefficients in forced vibration problems with the aid of Laplace transform, Computers and Structures 5 (1975) 263-269. G.V. Narayanan and D.E. Beskos, Numerical operational methods for linear dynamic systems. To appear. A. Papoulis, A new method of inversion of the Laplace transform, Q. Appl. Math. 14 (1957) 405-414. F. Durbin, Numerical inversion of Laplace transforms: an efficient improvement to Dubner and Abate’s method, Comp. J. 17 (1974) 371-376. T.V. K&m&n and M.A. Biot, Mathematical methods in engineering (McGraw-Hill, New York, 1940). S.H. Crandall, N.C. Dahl and T.J. Lardner, An introduction to the mechanics of solids, 2nd ed. (McGraw-Hill, New York, 1972). B.A. Boley and J.H. Weiner, Theory of thermal stresses (Wiley, New York, 1960). T.A. Cruse, A direct formulation and numerical solution of the general transient elastodynamic problem II, J. Math. Anal. Appl. 22 (1968) 341-355. R. Dubner and J. Abate, Numerical inversion of Laplace transforms by relating them to the finite Fourier cosine transform, J. ACM 15 (1968) 115-123. J.W. Cooley and J.W. Tukey, An algorithm for the machine calculation of complex Fourier series, Math. Comp. 19 (1965) 297-301. A.M. Ebner and D.P. Billington, Steady state vibration of damped Timoshenko beams, J. Struct. Div. ASCE 94 (1968) 737-760. G. WarzCe, Finite element method and Laplace transform -comparative solutions of transient heat conduction problems, Computers and Structures 4 (1974) 979-991.
MECHANICS COMPANY
AND ENGINEERING
21 (1980) 337-35.5
THERMALLY INDUCED VIBRATIONS OF BEAM STRUCTURES G.D. MANOLIS* Department
and D.E. BESKOS**
of Civil and Mineral Engineering, University Minnesota 55455, USA Received Revised
19 January 9 March
of Minnesota,
Minneapolis,
1979 1979
A general numerical method is developed for determining the dynamic response of beam structures to rapidly applied thermal loads. The method consists of formulating and solving the dynamic problem iq the Laplace transform domain with the aid of dynamic stiffness influence coefficients defined for a beam element in that domain and of obtaining the response by a numerical inversion of the transformed solution. Thus, the solution of the associated heat conduction problem, usually obtained by Laplace transform and needed for the computation of the thermal load, can be used in its transformed form. The effects of damping and of axial compressive forces on the structural response are also studied. Three examples are presented in detail to illustrate the proposed method and demonstrate its advantages.
1. Introduction The conventional thermal stress analysis of structures ignores the effect of inertia and, by using a predetermined time-dependent temperature distribution, essentially deals with a static problem where time is reduced to the role of a parameter. For slender structures under rapid surface heating however, such as beams, plates and shells, the inertia effects become significant and lead to thermally induced structural vibrations. Boley [l], [2] was the first to take structural inertia into account and analytically treat the problems of thermally induced vibrations of simply supported elastic beams and plates. His work was followed by the analytical studies of Kraus [3] on thermally induced vibrations of thin elastic spherical shells and of Yu [4] and Frisch [5] on observed thermally induced vibrations of flexible satellite booms subjected to variable heating. Boley [6], reconsidering the problem, derived a simple approximate formula for the ratio of the maximum dynamic to static deflection of heated beams and plates, and he then generalized it to include the effects of damping and of axial or in-plane loads. In addition to the aforementioned thermally induced vibrations in aerospace structures due to surface heating, similar vibrations due to radiation heating have been observed and studied in nuclear structures [7], [8]. A review article by Bargmann [9] discusses both types of thermally induced vibrations. This paper is concerned with thermally induced vibrations of structures composed of beams which are exposed to rapid surface heating. For complex beam structures, such as continuous beams and gridworks or frameworks, the analytical approach [l]-[5] based on the MindlinGoodman method [lo] for solving equations with time-dependent boundary conditions is not *Graduate Student. **Associate
Professor.
338
G.D. Manolis, D.E. Beskos, ~e~u~ly
induced oibrutio~ of beam structures
efficient, and numerical procedures should be used. In this work a general numerical method of analysis is developed which is capable of determining the elastic dynamic response of any plane beam structure to any time-dependent thermal loading. The method consists of applying the Laplace transform with respect to time to the governing equations of bending motion of an elastic beam element, solving for the transformed displacement function and constructing stiffness coefficients in the transformed domain. The problem is then formulated and solved in the transformed domain, and the transformed solution is inverted numerically to yield the dynamic response. The transformed thermal loading is usually obtained by solving the corresponding heat conduction problem in the Laplace transform domain. This method of solution of dynamic problems was first presented by Beskos and Boley [ll] and more recently generalized and extensively treated by Narayanan and Beskos [12]. The most important problem in this method is the numerical inversion of the Laplace transform, which is discussed in detail [12]. Here the method based on the trigonometric expansion as described by Papoulis [13] and the method based on sine and cosine transforms as described by Durbin El43 are utilized, The second method is more time-consuming than the first one but is more accurate for short, as well as long, time solutions. Other effects such as damping, rotatory inertia and shear deformation, as well as the effect of axial loading on flexure, can be easily considered by constr.ucting appropriate stiffness coefficients in the Laplace transform domain to include them. Here the effects of damping and of axial loads on the response are studied. Illustrative structural examples of beams and gridworks are presented to demonstrate the use and importance of the method, which presents distinct advantages over other well-known methods as far as accuracy, simplicity and efficiency are concerned.
2. General formulation and solution This section describes in detail the proposed numerical method for formulating and solving the general problem of thermally induced vibrations of beam structures exposed to rapid surface heating. Consider a uniform linear elastic Bernoulli-Euler beam element l-2 of length L, bending rigidity EI, cross-sectional area A and mass per unit of length m under rapid surface heating as shown in fig. 1. On the assumption that the resulting thermal loading is concentrated at nodes 1 and 2, the equation of flexural motion for the beam element l-2 is given by
where z, = u(x, t) stands for the lateral beam deflection along the positive y direction of fig. 1, and t denotes time.
Fig. 1. The flexural
beam element
in the Laplace
transform
domain.
G.L>. ~anol~s, D.E. Beskos, penalty
By following [Ill the Laplace transform initial conditions to yield
undyed vibrationsof beam secures
with respect to time is applied to (1) under zero
$+4zc4u‘=o, where a(x, S) =
I0
339
(2)
m U(X,t) e-“’ di
(3)
and 4K4 = ms2/EI. The
(4)
solution of eq. (2) is given in [15] as 6(x, S) = eKX[A COSKX) + B sin(&)] + eeKX[C cos(Kx) + D sin@)],
(5)
where K is given in terms of s by (4), and A, B, C and D are constants of integration. The beam element 1-2 has two degrees of freedom per node with positive directions for them as well as for their co~esponding nodal forces as shown in fig. 1. Thus, using 6(x, s) as the displacement function in a stiffness formulation, the beam element stiffness matrix can be constructed in the Laplace transform domain following standard procedures. The following relations are utilized in constructing the element stiffness matrix: 8(x, S) = -d;iildx, V(X, S) = -EI(d3E/dx3),
(6)
ti(x, s) = EI(d2$dx2), where 8 is the rotation, and v and $i are the shear force and the bending moment, respectively, with positive directions taken as in f16]. Finally the dynamic stiffness equation for the beam element l-2 in the Laplace transform domain takes the following static-like form:
(7)
In eq. (7), where the subscript T indicates thermal load, the dynamic stiffness influence coefficients in the Laplace transform domain are given by the following expressions: Dl, = 2NK2(y2 + 4ysc - l), 62, = Mq-y* d,,
+ 2y(l-
2s2) - 11,
= (4IvK”/vg[-y’(s3
lY,l = (2NK/%5)(s3
+ SC2+ c) - y(s3 + SC2- c)],
+ SC2+ s)(y - y2),
L522= N(y2 - 4ysc - l), n+* = (2rJ;tv7)[y’(s’
+ SC2- c) + y(s3 + SC2+ c)],
Is,, = -&,
D33 = Ds,,,
D43 = -LTZ,,
6,
= D2*,
(8)
340
C.D. Man~~is, D.E. Beskos, ~er~a~~y
induced vibrations of beam structures
where N = 2J51K/[y2 - 2y(1+ 2s2)+ 11, c = cos(KL),
s = sin(KL),
y = exp(2KL).
(9
With the introduction of the coefficients rS, the dynamic problem of a plane beam structure in flexural motion due to thermal loading is expressed in the Laplace domain in the static-like form F(S)T = ~(~)~(~),
(10)
where F(s)T and G(S) are the Laplace transformed thermal load and displacement vectors, respectively, while B(s) is the transformed dynamic stiffness matrix of the structure obtained by an appropriate superposition of element stiffness matrices as the one in (7). After application of the transformed boundary conditions, the transformed solution G(S) of eq. (10) is obtained, for a sequence of values of S, by a standard matrix inversion procedure usually performed numerically. From the above solution the response u(t) is obtained, for a sequence of values of t, by a numerical inversion of the Laplace transformed displacements. Heating of a beam element produces a thermal nodal loading consisting in general of thermal bending moments and axial forces. Thermal axial forces can be disregarded on the assumption that axial deformations are not restrained, while thermal moments are given by MT=Ea
Ty W
(11)
A
where (Yis the coefficient of thermal expansion, which is assumed to be constant, and T is the temperature dist~bution in the beam, in general a function of X, y and t. The temperature dist~bution T is obtained by solving the heat conduction problem defined for the heated structure. Consider an isotropic solid body with a constant thermal diffusivity U, mass density p and specific heat c, subjected to arbitrary thermal boundary and initial conditions on its surface and with internal heat generation 4 per unit time per unit volume, which is of importance in problems of radiation heating [7], [8] but is taken as zero for problems of surface heating as in the present case. The general equation of heat conduction is (see [17])
Uv*T+-g, PC
(12)
where T is the temperature of the body, varying with space and time, and V2 is the Laplace’s operator_ Eq. (12) is supplemented with appropriate boundary and initial conditions to form a well-posed boundary value problem. The five principal boundary conditions of heat conduction are described in detail in [17]. For the problem of surface heating of a uniform flexural beam element of rectangular cross-section with depth h and width b one has that 4 = 0 and T = T(y, t) in (12) so that this equation becomes (13)
G.D. ~a~olis,
D.E. Beskos, verbally
induced uibraiio~s of beam sfruc~ures
341
If 7’ is a function of x as well, one can afways subdivide the structure into elements with constant temperature along the beam element length. In the present problem the initial condition is assumed to be
0) = 0,
VY,
(14)
while the boundary conditions K
are
5 =-Q(t)
at
KO
y = h/2,
at y =
ay-
(15)
-h/2,
indicating application of a heat input Q(t) on the top beam fibers and insulation of the bottom ones. The constant K in (15) is the thermal conductivity of the material and is related to u by
One usual way of solving the heat conduction boundary value problem consisting of (13)-(H) is by the method of Laplace transform. With the application of the Laplace transform with respect to time eqs. (13~(15) are transformed into the system (d2ndy2) - (s/u)l?; = 0, K(dmdy) = -b(S),
at
y = h/2,
dp,‘dy = 0
at
y = -h/2,
F = F(y,
s) = fin Tfy,
t) e-“’ dt,
&)
= lrn Q(t) e-“’ dt,
(17)
0
which yields the following solution for the temperature F(y, s) = A eAy+ B eeAy,
in the Laplace transform domain: (18)
where h = %5/U, A = -(~/~~)~e~~~~/(lB = Since the transform by Laplace of bending
-(@hK)[ehh’2/(1
-
eU”)],
(19)
e2Ah)].
problem of the structural vibration is formulated and solved in the Laplace domain, there is no need to obtain the temperature distribution in the time domain inversion. In view of (ll), the nodal thermal loading in the Laplace domain consists moments which for a rectangular beam are of the form M;_(s) = =b
1-y; ?‘(Y, 4~ dy,
where F(y, s) is given by (18) and (19).
(20)
342
G.D. Manolis,
D.E. Beskos,
Thermally
induced vibrations of beam structures
The heat input Q(t) can be any fast-rising function of time such as or a gradually applied (ramp-type) one with small rising time. In nonlinear function of time which can be approximated by a piecewise case, of course, it is preferable to compute o(s) and consequently following [ 121.
a suddenly applied one fact, Q(t) can be any linear function. In that MT(s) numerically by
3. Numerical inversion of Laplace transform Let f(t) be a real function of t, with f(t) = 0 for t -=L 0. The Laplace transform f(s) of f(t) is defined by
T(s)= 6 f(t)
e-“’ dt,
(21)
and its inversion formula is given by
where p > 0 is arbitrary but is greater than the real part of any of the singularities of f(s), and s is a complex number with Re(s) I j3 > 0. In many engineering applications f(s) is either too complicated to be inverted analytically by evaluating the integral of (22) or is given in numerical form. In those cases a numerical inversion of the Laplace transform is imperative. Reference [12] discusses this problem in detail and investigates the accuracy, computational efficiency and range of applicability of eight methods of numerical Laplace inversion as applied to dynamic problems. Among those methods two are briefly presented and used in this work. The first method is due to Papoulis [13] and determines f(t) in terms of the values of f(s) on a finite sequence of equidistant points s, on the real axis by the formula f(o) = 5 C, sin(2n + l)O, n=O
(23)
where N is a finite positive integer number, the coefficients C,, are given in [13] as functions of the values of T(s)), and 0 = cos-‘(e-“‘),
(24)
where S > 0 is a constant selected by the user. This method has the advantage that it works with real data and provides easy-to-use expressions for the evaluation of the C”, thus requiring small amounts of computer time. Unfortunately, the accuracy of the method is only acceptable for early times, and improvement of its accuracy by increasing hi is not possible because the computation of the C, becomes an unstable process for large N [18]. The second method of inversion is due to Durbin [14] and is actually an efficient improvement
G.D. Manofis, D.E. Beskos, Thermally induced vibrations of beam structures
343
of the method of Dubner and Abate [19], which is based on the finite Fourier cosine transform. Durbin [14] combined both finite Fourier cosine and sine transforms to obtain the inversion formula f( tj )
x
eaiAr [-+Rev@)]+Re{y
(A(n)+iR(n))W”]], tl=O
where
Re{f@
+ i(n + 1N) F
B(n) = 5 Im{f(P I=0
+t i(n + fN)+
A(n) = i
1, tz = 0, 1, . . . )
and f(t) is computed for N equidistant points jAt = jT/N, j = 0, 1,. . . , N - 1. It is suggested that for L x N ranging from 50 to 5000 one should select fit = 5 to 10 for good results, where t^ is the total time interval of interest. The computations involved in (25) are performed by employing the fast Fourier transform algorithm of Cooley and Tukey [20]. The accuracy of Durbin’s method [14] is very high even for late times, and its only disadvantage is the fact that it works with complex data and thus is more time-consuming than the method of Papoulis [13]. This disadvantage is however offset by its high accuracy, as the numerical examples in this work will clearly demonstrate.
4. Damping and axial force effects The effect of damping on the thermally induced motion of a beam has been investigated by Yu f4] and Boley [6]. Yu [4] considered internal viscoelastic damping of the beam material by employing the constitutive law of the standard linear solid and installing a viscous-fluid damper at the tip of the cantilever beam when modeling a satellite boom. Boley [6] studied the damping effect by considering viscous damping in a single degree of freedom system approximating beams or plates. In this work both internal viscoelastic damping and external viscous damping are considered. Internal viscoelastic damping is accounted for by assuming that the beam material obeys the constitutive law of a Kelvin viscoelastic model: CT= E(E + j-k),
(27)
where u is the stress, E is the strain, g is the strain rate, E is the moduius of elasticity, and f is the internal damping coefficient of the material (which is usually very small for metals). Application of the Laplace transform as defined by (21) under zero initial conditions on (27) yields (T = E(l + fS)Z,
(28)
G.D. Manolis, D.E. Beskos. Thermally induced vibrations of beam structures
344
where Cr and E are the transformed stress and strain. Eq. (28) indicates that internal viscoelastic damping can be very easily taken into account if E is replaced by E(l + fs) in eqs. (9) and (20). In the case where there is external viscous damping eq. (1) is replaced by (29) where c is the viscous damping coefficient. Application defined by (3) yields
of the Laplace transform
to (29) as
$$+4K4C=0,
(30)
4K4 = (ms’+ cs)/EI.
(31)
where
Thus, in view of the fact that (2) and (30) are identical in form, external viscous damping can also be very easily taken into account if (4) is replaced by (31). The effect of axial loads or beam-column effects on the thermally induced motion of a beam has been approximately treated by Boley [6], who utilized the fact that, for simply supported beams under certain types of loads, the magnification of the static or dynamic deflection under an axial compressive load P can be approximated by the formula 1
& =1 - (PIP,,)’
(32)
where PC, is the buckling load. In this work the effect of axial loads on the response of beam structures to thermal loading is determined in an exact way by constructing the coefficients Dii which include the effect of axial loads on flexure. In the presence of a constant compressive force P, eq. (1) is replaced by EI ~+P$+m$=O, a4v which, after application
(33)
of the Laplace transform as defined by (3) becomes (34)
where 2~ = PIEI,
v = ms2/EI.
(35)
The solution of (34) is given by a(~,
s) = A e*‘”+ B e’2x+ C eAp + D e”“,
(36)
G.D. Manolis, D.E. Beskos, Thermally induced vibrationsof beam structures
where A, I?, C, D are constants numbers given as
of integration,
345
and AI, AZ, hJ, A4 are in general complex
A,=&-,
A2 = -J-CL
+ e,
A3=+p-m,
A4 = -d-p
-m.
(37)
On the basis of the displacement function (36) the D, including the effect of P can be constructed as in section 2. Here, however, due to the lengthy algebraic computations involved, explicit expressions for the Dij are not derived. Only the constants of integration A, B, C and D of (36) are explicitly determined and given in the appendix for the necessary displacement states which define the Dij. For a beam-column element the Dij are computed numerically on the computer with the aid of (36) and v(x, s) = -EI(d3G/dx3) - P(dG/dx), n;i(x, s) = EI(d2G/dx2).
(38)
The case of a tensile axial force P can be very easily taken into account by simply replacing P by -P in all the previous expressions. The combined damping and beam-column effects can also be treated on the basis of eqs. (33)-(38) either by replacing E by E(l + fs) for the case of internal viscoelastic damping or by defining v as v =
(ms2+ cs)/EI
(39)
for the case of external viscous damping.
5. Structural examples This section describes the detailed solution of three structural examples in order to illustrate the use of the proposed numerical method, demonstrate its advantages and provide qualitative as well as quantitative results for the problem of thermally induced vibrations of plane beam structures. Example 1 Consider a uniform simply supported beam of length L, bending rigidity EI, mass per unit of length m and of rectangular cross-section with width b and depth h as shown in fig. 2. A
Fig. 2. A simply supported
beam under suddenly
applied
heat input.
346
G.D. Manolis,
D.E. Beskos.
Thermally
induced
vibrations of beam structures
step heat input Q(t) of intensity QO, uniformly distributed along the length of the beam, is assumed to be applied along the edge y = h/2, while the edge y = -h/2 is insulated, as fig. 2 shows. The deflection and bending moment of this beam have been analytically obtained by Boley [l], who presented his results in terms of the following dimensionless parameters:
T =
utlh2,
VW
V(& T) = 7r4~v/(192Q~L2), mT(& T) =
(4Oc)
n4~hfT/(192EI@),
(404
where z, and MT stand for deflection and thermal moment, respectively, and all the other quantities which appear in (40) have been defined in section 2. It was shown in [l] that the governing nondimensionless parameter of the problem, which is an index of the importance of the inertia effect, is B. This is defined by (40a) as the square root of the ratio of the characteristic thermal time h2/u to the characteristic mechanical time (-the first natural period of vibration). It was found in [l] that the importance of the inertia effects increases with decreasing B so that for B = 0 they actually prevent any beam deflection. Also, as B approaches infinity, the inertia forces disappear and the static solution results. Furthermore, it was also shown in [l] by a quantitative analysis on an aluminum beam that small values of B for which inertia effects are important correspond to very thin beams. In this work, on the assumption that the beam material is aluminum, the following numerical data was adopted for the computation of the midspan deflection and bending moment of the beam l-3 of fig. 2: b =
1 in,
h =
0.1539 in
L = 40 in,
m = 4 x 10d5 lb sec2/in2,
E = 10.66 X lo6 lb/in2,
u = 0.1333 in21sec,
K
cy = 12.25 X 1O-6l/“F,
Q. = 1 Btu/in2 sec.
=
0.0027 Btu/in set “F, (41)
On the basis of the above data one has that B= 1,
T = 5.62789&
V = O.O699v,
mT = O.O34966M,
(42)
In order to determine the midspan dynamic deflection n2 of the beam l-3 of fig. 2, a node is introduced at midspan, thus dividing the beam into two elements. Following the general procedure of section 2 and utilizing the boundary conditions & = iT3= 0, 0, = 0 and 0, = -&, the deflection C2 due to nodal thermal loads M1, = -Z&, = i& and l&- = 0 takes the form
G.D. Manolis, D.E. Beskos, Thermally induced vibrations of beam structures
347
where the Oij are given by (8), and A& is given by (18)-(20) for a(s) = Qols. The response ~)~(t)was obtained by numerically inverting the Laplace transformed solution U2given by (43). The methods of Papoulis [13] and Durbin [14], as described in section 3, were both employed, and the results are shown in fig. 3 together with the exact analytical results of [l]. It is apparent in fig. 3 that, while Durbin’s results indicate a perfect agreement with the exact solution, this is not the case with Papoulis’ results, which present a reasonable agreement only for the first cycle of vibration. The CP computer time spent on the CDC Cyber 74 computer of the University of Minnesota Computer Center was 2.118 set for Durbin’s and 1.285 set for Papoulis’ methods, respectively, for 250 time points. The values of 6 = 5.866 and N = 9 for Papoulis’ method and L = 1, N = 250, t^= 1.0 set and pt^ = 6 for Durbin’s method were used in this example. Due to the high accuracy results of Durbin’s method all further numerical Laplace inversions were performed by that method, even though more computer time is required by that method _than by Papoulis’ method. The bending moment at the midspan of the beam l-3 of fig. 2 was also computed by observing that
M,=M;i-AT;,
(44)
where the dynamic moment A?,” is obtained by considering the beam element l-2 and solving for the internal element forces in terms of known nodal displacements and has the form
and the static moment A?; is equal to fi, where A& is given by (18)-(20) for a(s) = QJs. The bending moment A& given by (44) was inverted numerically by Durbin’s method, and the results show excellent agreement with the exact ones of [l], as fig. 4 demonstrates in a dimensionless plot of m2 = 7r4~M2/192QcuEI versus 7 = ut/h*.
-
Bolcy PI
--*-*
Papoulls
0000
Durbln
[Ia] [I41
Fig. 3. Dynamic midspan deflection of the simply supported beam for B = 1, f
= 0.
348
G.D. Manolis, D.E. Beskos, Thermally induced vibrations of beam structures
m
-
a3
B&y
* * 0 The
[I] method
02
Fig. 4. Dynamic midspan bending moment of the simply supported beam for I? = 1, ,f = 0.
The effect of internal viscoelastic damping on the midspan deflection is portrayed by fig. 5. The computations were done on the basis of (43), with E being replaced by E(l +fs) everywhere. Because the internal damping coefficient f is very small for metals and ranges between 0.005 and 0.020 (see [21]), the values 0.005, 0.010 and 0.015 were used in this example. Fig. 5 shows clearly the attenuating effect of damping on the dynamic response, which increases rather rapidly with increasing values of f and time T. Eventually, after some time the damped dynamic deflection approaches the static one which corresponds to the case of zero inertia effects. This indicates that damping is in general an effective way to reduce inertia effects in a thermal stress structural problem. 4 V
f
0.06 -
0
0.5
1.0
1.5
-c
Fig. 5. Dynamic midspan deflection of the simply supported beam for various values of the intertia index B and the damping f.
G.D. Manolis, D.E. Beskos, Thermally induced vibrations of beam structures
349
A
v,
(If-l)
3.0 -
This
method
2.0 -
1.0 -
c
0
Fig. 6. Dynamic
midspan
deflection
0.05
0.20 t
0.15
0.10
of the simply supported
(Set)
beam for B = 1, f = 0 and various
values of P/P,,.
Finally fig. 6 depicts the amplification of the midspan deflection v2 of the beam of this example due to the presence of an axial compressive force P
+
t
Q(t) QO
Q(t)=0
\/
1
i/l
1 4 l_ A
El m T
+-C--+-C_, 2 Fig. 7. A two-span
continuous
El m y 2
beam under suddenly
applied
heat input.
G.D. Manolis, D.E. Beskos, Thermally induced vibrations of beam structures
350
obtain the following expressions for the rotations at the supports of the continuous
beam:
-0, = -Ii&& &=o,
(46)
where fir is given by (18j(20) for o(s) = Qofs, and & is given by (8). A numerical Laplace inversion of 6, yields o,(t), which is plotted in fig. 8 for f = 0 and f = 0.015. Fig. 8 shows that the inertia effect is not so significant in this case as B,(t) oscillates with small amplitude about its static value. In fact, for a damping value of f = 0.015 there is virtually no oscillation, and o,(t) rather rapidly approaches the static solution.
Fig. 8. Dynamic rotation @x(t)of the continuous beam for B = 1, f = 0 and 1.5%.
Example 3 Consider a plane horizontal square gridwork simply supported at its four corners as shown in fig. 9. The gridwork consists of beam elements of length L, = 24 in and depth h = 0.2 in and with all other geometrical and material properties as in (41). It is assumed that a suddenly applied heat input acts on the top fiber of all the beams of the gridwork, while their bottom fibers are insulated. It is further assumed that axial and torsional deformation of the beams can be neglected so that element equilibrium equations in the Laplace transform domain are of the form (7). In this example, due to the large number (23) of free degrees of freedom of the gridwork, the formulation and solution of the problem was done exclusively by computer. A general computer program was written in FORTRAN capable of performing a complete dynamic analysis of a gridwork under thermal loading on the basis of the proposed method as described in section 2. The program constructs the total dynamic structural stiffness matrix fi and solves equation (10) numerically by Gauss-Jordan elimination for a sequence of complex values of s. The
G.D. Manolis, D.E. Beskos, Thermally induced vibrationsof beam structures
351
Fig. 9. A plane horizontal square gridwork under suddenly applied heat input.
transformed solution thus obtained in discrete form is inverted numerically by Durbin’s method to yield the dynamic response. Using this program, the vertical deflection at node 9 of the gridwork of fig. 9 was determined as a function of time, and the result is shown in fig. 10. The total CP computer time for this computation was 20 sec.
20-
1.0-
i
00
0.1
0.2
0.3
t mx)
Fig. 10. Dynamic vertical central deflection of the gridwork for B = 1, f = 0.
6. Concl~ions On the basis of the preceding discussion one can draw the following conclusions: 6.1
A general numerical method for determining the thermally induced dynamic response of elastic beam structures has been developed. The method formulates the dynamic problem in the Laplace transform domain-in a static-like form by using Laplace transformed dynamic stiffness influence coefficients -ct, and solves it (there) numerically. The dynamic response is finally obtained by a numerical inversion of the transformed solution.
352
G.D. Manolis, D.E. Beskos, Thermally induced vibrations of beam structures
6..2 Only Laplace transformed dynamic thermal loads are used in this method. Thus the associated transient heat conduction problem under any kind of boundary conditions has to be formulated and solved in the Laplace domain only, and there is no need for its time domain solution. 6.3 In the framework of the linear theories employed here the use of the Dij coefficients leads to the “exact” solution of the dynamic problem since the displacement function used for their construction is the exact solution of the transformed equation of motion. Thus the present method can provide a basis for comparing the accuracy of other approximate methods such as the finite element method or the method of Boley [6], which is based on the assumption that the first-mode response is legitimate for many practical problems. 6.4 Various effects on the response of beam structures such as damping, axial forces, axial or torsional deformations, rotatory inertia and shear deformation can be very easily taken into account by constructing appropriate Dij coefficients which include these effects. In this work the effects of damping and axial forces were considered and studied. As far as damping is concerned, use of Dij coefficients permits consideration of different amounts of damping in different beam elements, thus enabling one to effectively control the response by appropriate damping changes in particular structural elements. 6.5 The proposed method is quite general and capable of solving any dynamic problem of beam structures and not just thermally induced vibrations. From this point of view it can be considered as a generalization of the method proposed by Beskos and Boley [ll]. Because in problems of thermally induced vibrations the principal results for beams and plates are almost identical regardless of the aspect ratio of the plate or the rapidity of the heat application and even of the support conditions [6], the results obtained by the proposed method for beam structures can be quite informative about the response of their associated plate structures. Apart from that, the basic ideas of the method can be extended to treat dynamic problems of plate and shell structures directly. 6.6 The proposed method appears to be better than the conventional finite element method in conjunction with either modal analysis or numerical integration because i) it provides the exact solution of the dynamic problem and not an approximate one as in the finite element method (which discretizes the continuously distributed mass of physical beam elements), and ii) it considerably simplifies the solution of the associated heat conduction problem. Furthermore, the method does not require prior knowledge of natural frequencies and modal shapes as in modal analysis. The only disadvantage of the method has to do with the numerical Laplace
G.D. Manolis, D.E. Beskos, Thermally induced vibrations
of beam
structures
353
inversion, where one has to decide between relatively low accuracy with less computer time and high accuracy with more computer time and select his inversion algorithm on this basis. Reference [12] indicates that, for a given dynamic input, numerical integration techniques such as Newmark’s p-method or Wilson’s O-method require in general less computer time for a dynamic problem than numerical Laplace transform techniques for the same degree of accuracy. The latter techniques however have the characteristic properties that they permit one to perform numerical inversion for just one nodal displacement of interest independently of the others, or determine a particular nodal displacement at a desired time without prior knowledge of its values at all previous times. From these view points numerical Laplace transform techniques are more efficient than numerical integration methods. Also, as [22] points out, Laplace transform techniques based on algorithms of numerical inversion working with real data seem to be more economical than numerical integration for large systems to be studied over a large time domain because no more than 20 system resolutions need to be computed.
Acknowledgement
The authors are grateful to the University of Minnesota Computer Center for making its facilities available to them and to Dr. F. Durbin for supplying them with his numerical Laplace inversion computer program.
Appendix. Determination displacement states
of the constants of integration A, B, C and D of (36) for the various
1st displacement state. 6(0, s) = 1, ~‘(0, S) = I?(& s) = IT’& S) = 0, where the prime indicates differentiation with respect to x. For this case the use of (36) results in the system A+B+C+D=l, hlA + A,B + h,C + AbD = 0, ehlLA + eAzLB + eA*C + eALD = 0, Al
eAILA + A2 ehzLB + A3 eA*C + A, e”&D = 0,
(A4
which has the solution A = (l/A)[hz(A4 - A3)e(*3+*4)L + A3(A2- A4)e(h2+A4)L + A4(A3- AZ)e(“2+A3’L], B = (l/A)[A1(A3 - A4)e(h3+A4)L + A3(A4- A,) e(h1CA4K + A4(Al- AS)eQ’+*3)L], C = (l/A)[AI(A, - AZ)e(A2+h4)L + A2(Al- A4)e(h1+*4)L + A4(Az- AI) e(A*+*2)L], D = (l/A)[A,(A,-
As) e(A2+*3w + A,(A, - A,) e(A1+*3)L + Ag(AI- AZ)e(h’+*2)L],
64.2)
where A = sAlA + (A, - A3)(A4- AZ)[e(AlfA3ZL + e(A2+r4)L] + (A, - A4)(A2- A3)[e(A1+h42L + e(A2+*3)L].
64.3)
354
ti.D. Manolis, D.E. Beskos, Thermally induced vibrations of beam structures
2nd displacement state. a(O, s) = 0, U’(0, s) = 1, U(L, S) = c’(L, S) = 0. For this case the use of (36) results in the system A+B+C+D=O, A,A + h2B + h,C+ h,D = 1, eAILA + eAzLB+ eAXLC+ e*&D = 0, A,
e”lLA + A2 eA2LB + A3eA’LC +
Aq e**D
(A.4
= 0,
which has the solution
A =
(l/A)[(A,
-
B
=
(l/A)[(A, - A3)e(*3+*4jL + (A, - A4)e(*l+*d)L + (A3_ A,) e(*l+~T)~],
C
=
(l/d)[(hz-
D
=
(l/A)[(A, - A*)e(A*+AX)L + (A, - AX)e(*l+*J)L+ (A, _ A,) e(*~+~~)~],
Ad)
e(A3+A4)L + (Ad-- hZ)e(*2+*4)L
+ (AZ-
~~)~f*z+AdL],
Aq)e(A2+A4)L + (A4- A,) e(A1+*4)L + (A, _ AZ)e(*~+~~)~~,
(A.3
where A is again given by (A.3). Eqs. (38) computed at nodes 1 and 2 of the beam element of fig. 1 become V(O,s)=-EI(AA:+BA;+CA:+DA:)-P(AA1+BA,+CA,+DA,), M(o,s)=EI(AA:+BA;+~A:+DA:), v(L, s) = -EI(AA
f e*IL +BA~eA2L+CA~eA~+DA~eA4L)
- P(AA, eAIL+ BAZeAzL+ CA, e*& + DA, e*&), I@(& s) = EI(AA: eAIL+ BA$ eA2L+ CA: e*$ + DA: e*&).
(A.6)
Using (A.6), the Oij coefficients for the 1st and 2nd displacement states can be obtained by utilizing the values of the constants A, B, C and D for the corresponding displacement state as given by (A.2) and (AS). Thus the 1st ,and 2nd displacement states in conjunction with the proper sign convention consideration provide the first and second columns of the beam element stiffness matrix 0. The third and fourth columns can be very easily obtained from the first two by using the Betti-Maxwell reciprocal relations. References [l] [2] [3] [4]
B.A. Boley, Thermally induced vibrations of beams, J. Aerospace Scs. 23 (1956) 179-181. B.A. Boley and A.D. Barber, Dynamic response of beams and plates, J. Appl. Mech. 24 (1957) 41Wl16. H. Kraus, Thermally induced vibrations of thin nonshallow spherical shells, AIAA J. 4 (1966) 500-505. Y.Y. Yu, Thermally induced vibration and flutter of a flexible boom, J. Spacecraft and Rockets 6 (1969) 902-910. [S] H.P. Frisch, Thermally induced vibrations of long thin-walled cylinders of open section, J. Spacecraft and Rockets 7 (1970) 897-905.
G.D. Manolis, D.E. Beskos, Thermally induced vibrations of beam structures [6] B.A. Boley, [7] [8] [9] [lo] [ll] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22]
355
Approximate analyses of thermally induced vibrations of beams and plates, J. Appl. Mech. 39 (1972) 212-216. W.C. Lyons, Comments on heat induced vibrations of elastic beams, plates, and shells, AIAA J. 4 (1966) 1502-1503. H. Bargmann, Dynamic thermal shock resistance, pp. 174-181, in: J.L. Zeman and F. Ziegler (Eds.), Topics in applied continuum mechanics (Springer, Berlin, 1974). H. Bargmann, Recent developments in the field of thermally induced waves and vibrations, Nut. Eng. Design 27 (1974) 372-385. R.D. Mindlin and L.E. Goodman, Beam vibrations with time-dependent boundary conditions, J. Appl. Mech. 17 (1950) 377-380. D.E. Beskos and B.A. Boley, Use of dynamic influence coefficients in forced vibration problems with the aid of Laplace transform, Computers and Structures 5 (1975) 263-269. G.V. Narayanan and D.E. Beskos, Numerical operational methods for linear dynamic systems. To appear. A. Papoulis, A new method of inversion of the Laplace transform, Q. Appl. Math. 14 (1957) 405-414. F. Durbin, Numerical inversion of Laplace transforms: an efficient improvement to Dubner and Abate’s method, Comp. J. 17 (1974) 371-376. T.V. K&m&n and M.A. Biot, Mathematical methods in engineering (McGraw-Hill, New York, 1940). S.H. Crandall, N.C. Dahl and T.J. Lardner, An introduction to the mechanics of solids, 2nd ed. (McGraw-Hill, New York, 1972). B.A. Boley and J.H. Weiner, Theory of thermal stresses (Wiley, New York, 1960). T.A. Cruse, A direct formulation and numerical solution of the general transient elastodynamic problem II, J. Math. Anal. Appl. 22 (1968) 341-355. R. Dubner and J. Abate, Numerical inversion of Laplace transforms by relating them to the finite Fourier cosine transform, J. ACM 15 (1968) 115-123. J.W. Cooley and J.W. Tukey, An algorithm for the machine calculation of complex Fourier series, Math. Comp. 19 (1965) 297-301. A.M. Ebner and D.P. Billington, Steady state vibration of damped Timoshenko beams, J. Struct. Div. ASCE 94 (1968) 737-760. G. WarzCe, Finite element method and Laplace transform -comparative solutions of transient heat conduction problems, Computers and Structures 4 (1974) 979-991.