SOLUTION OF TIME DEPENDENT BOUNDARY VALUE PROBLEMS BY THE BOUNDARY OPERATOR METHOD H. D. FISHER, M. M. CEPKAUSKASand S. CHANDRA Reactor Engineering Department, Combustion Engineering Incorporated, Windsor, CT 06tl95,U.S.A. (Received13 September1978;in revisedform 14 December1978) Abstract-A new procedure, the Boundary Operator Method (BOM), for deriving normal mode response formulas in systems with time dependent boundary conditions is demonstrated. The method of analysis is applied to a f&e, nonhomogeneous bar. This example illustrates the following advantages of the BOM: (1) The resulting formulas for the displacement and stress are valid for the entire history of the system and for Dirichlet, Neumann or Cauchy boundary conditions. (2) The uniqueness of the solution is unobscured since, unlike the Mindlit&oodman method, the present solution technique does not require spe$fic spatial functions which render the transformed boundary conditions homogeneous. (3) The evaluation of the solution is expedited by the absence of integrals containing auxiliary spatial functions. (4) The analysis employs standard mathematical techniques.
NOMENCLATURE At, 4, Ck constants, F”L”p CO = (&J/P) elastic wave velocity, FOL’T-’ cm constant, F%?P Di cos ai + sin qi(d/&), FL”‘P fi dimensionless temporal function, F”Lo’P gi dimensionless spatial function, F’L’? I, = ” XL2de* FLOP I l, L length of bar, FL?” K, upper limit of summation in eqn (29), F”Lop real exponent boundary operator detined in eqn (l3), FL”T” Tk dimensionless temporal function, FL”T” real time FOL’TO U displacement in X direction, FLIT’ xk dimensionless eigenfunction, F”Lop X axial coordinate, F’L’T’ d$t$onless real numbers, FL”7’a
ol t
FLOT“ dimensionless axial coordinate, FL”T” dimensionless boundaries of bar, FL’T” Cl, c2 dimensionless displacement, FL”T” I K dimensionless coefficient F”Lop eigenvabres, FLT ’ Ak P mass density, F’LdT2 dimensionless amplitude of stress pulse,F”LoT+’ UO 7 = C&L dimensionless time F”LoTo _ -4 = UIKL dimensionless disoiacement. F”Lop * denotes dummy variable I denotes derivative wrt c denotes derivative wrt 7 C= HKL
INTRODUCTION
variety of computational techniques have been employed to analyze one dimensional, finite, continuous systems subjected to time dependent boundary conditions. Cinelli[l] has ascertained the transient response of thick elastic cylinders and spheres to dynamic surface loadings by employing a finite Hankel transform. After a transformation that resulted in homogeneous boundary conditions, Lee[2] also invoked the finite Hankel transform to study wave propagation in a nonhomogeneous rod. Schreyer [3] has developed an inverse procedure that yields exact solutions to the one-dimensional inhomogeneous wave equation. Reference [3] employed A
60?
608
H. D. FISHER, M. M. CEPKAUSKAS and S. CHANDRA
these transformation relations to compute amplitude amplification vs dimensionless frequency at various locations in a bounded gradient layer resting on sinusoidally excited bedrock. The method of characteristics was applied to the class of problems under discussion in Chou and Greif[4] and Greif and Chou[5]. In AIzheimer and Forrestal[6], the circumferential stress history at the inner surface of several thick cylinders subjected to a step pressure was calculated using the Laplace transform technique in conjunction with asymptotic expansion of the transform variable. Good agreement for the early time response was demonstrated between these results and the numerical studies of [4]. Lindholm and Doshi[7] utilized the method of virtual work to quantify the response of an elastic, nonhomogeneous, finite bar to a pressure pulse. This procedure was also employed by Fisher[8] in a study of wave propagation in a hollow conical frustum. Recently, Pao and Ceranoglu [9] have studied thick-wall spherical shells using ray theory. Numerical results presented in [9] quantify the amplification and the compression-tension oscillation of the radial stress history at various points in the hollow sphere. A general method that’ is applicable to all the problems discussed above is described in Mindlin and Goodman[ lo]. In addition to generality, the algorithm contained in [lo] has the advantages of being readily understood and yielding solutions valid throughout the time of interest, i.e. for temporal values which preclude asymptotic expansion of the Laplace transform. These qualities have resulted in wide application of the Mindlin-Goodman method as indicated in Meirovitchllll, Epstein[I2] and the bibliography in [12]. The essence of the procedure is to extend the method of separation of variables to the solution of partial differential equations governing continuous systems with time dependent boundary conditions. This is accomplished in a system with one spatial variable by expressing the original dependent variable as the sum of a new dependent variable and a finite series. Each term in this series is the product of a boundary condition and an auxiliary function whose argument is the spatial variable, The only constraint on the auxiliary functions is that they contain a sufficient number of terms to render the transformed boundary conditions stationary. Due to the homogeneity of the resulting boundary conditions, the transformed system of equations is amenable to solution by separation of variables. The desired solution of the original time dependent boundary value problem is then obtained from the stated transformation equation. As inspection of the solution given in [lo] for a vibrating beam restrained by time dependent boundary conditions reveals, the spatial dependence of the auxiliary functions required in a given problem must be explicitly stated and the integrals containing these functions must be evaluated to obtain the desired solution. Reflection on the technique outlined above results in the conclusion that the introduction of auxiliary functions is a mathematical convenience, and that the possibility exists of deriving solutions that do not contain such functions. This follows since well posed problems in partial differential equations are known to have unique solutions, i.e. solutions that are independent of the particular auxiliary functions utilized. For one-dimensional wave propagation in a finite, nonhomogeneous bar with time dependent displacement and/or stress boundary conditions, this paper demonstrates a technique, the Boundary Operator Method (BOM), that eliminates the auxiliary functions from the solutions for the displacement and stress in the bar. The generality of the BOM is identical to [lo] and, for this reason, the BOM can also be applied to one-dimensional wave propagation problems involving beams, spheres and cylinders composed of homogeneous or nonhomogeneous material.
ANALYSIS
Based on [7], the equation of motion for the free longitudinal motion in a thin bar with a lengthwise variation in elastic modulus is:
(1) where WC T) is the dimensionless displacement from the initial position, n is an arbitrary power that describes the axial modulus variation and e and T are the dimensionless spatial and temporal variables, respectively (see Nomenclature). Generalized time dependent boundary
Solution of time dependent boundary value problems by the boundary operator method
609
conditions are described by: Di[Y(eiy T)] = Y(Ci, 7) COSO!i + $(ei,?)
sin Qi =fi(T); i = 1,2
(2)
where ~1 and ez correspond to the left and right boundaries of the bar, respectively, and the ai are real numbers. The associated initial conditions are:
$
Y(c, 0) = Ye(E),
(E,0) = *o(E).
(3)
$,gi(E)fi(T)-
(4)
Now define the transformation:
*(G T)= [(‘S’T)+ With the aid of (4), eqns (l)-(3) are written as:
Di[l(Ei, T)] =
f,(T)-,$ Di[gj(Ei)lfj(T)ii = 132
(6)
Equation (6) is made homogeneous by requiring that: Di[&(Q)l=
&j7
i=1,,2;j=l,2
(8)
whereSii=Ofori#jandSri=lfori=j. A solution to (5) is now sought in the form:
(9) where the Xk(ci) Satisfy e”i%$+
ne”-‘xk = -Azxk,
Di[xk(6)1
= 0
(10)
and the rigid body eigenfunction, X0, has been omitted from the summation. Substituting (9) into (5), multiplying both sides by Xk, integrating from cl to er and utilizing the orthogonahty of the eigenfunctions yields:
(11)
(8) and (lo), the equations: gj(~i))Xi(~i)-gj(Ei)Xk(~i) = [-xk(ei) Sill (Yi +Xi(Ei)
COS
ai]&,, i = 1,2; j = l-2
(12)
610
H. D. hfE&
M. M. CEPKAUSKAS and S. t%ANDRA
are derived. Thus integrating the tist term in the numerator of (11) by parts, invoking (12) and defining the boundary operator:
Oi[Xk(Ei)]= (-l)i(Ei)“[Xk(ei)
Sill aj
O&%(cii)l= (- l)‘IX,(eJ sin ai -XL(q)COS
-X;(Q) ai],
1,2; Ei# 0 n = 0 and one of Ei= 0
COS (Yi}, i =
i =
1,2;
(13)
simplifies (11) to: (14) The solution of (14) is:
Integrating the last term
in
the
numerator
of (IS) by parts
(15) to be
permits
written
as:
sin A~(T - T*) dr*
T&(T)= AL cos ALT+ & sin At? f (A&)-’ ~‘~i~X~(~~~~lfi(r*)
\‘*g&k de*)
+ g(O) sin At.7 + A~i(O) COS ANT - Adi(
s (16)
Cf
Hence, . Tk(7)
=
-AkAkSill At7 + A&
AkT + 1k-I 2
COS
(oi[&(ei)]
i- Gi(O)COS From
(7)
[
fi(T*)
Ad:(O)Sin
hkT -
Ak(T- r*) dr*
COS
I\kT -f;(T)]
c&!&k
de*)*
(17)
and (9)
Specifying r = 0 in (16) and (17) and comparing with (18) yields: Ak = &
TttT)=zk-‘(I(
=
I
zk-*
(r\&)-*
de*
"[~~~*)-,~g,(~*)~(o)]xk (I.
Vlo(C*) COSAk* -br\k-“@o(6*) Ski hkT + Ak-*2
de*
&(~*)ji(o)]Xk
[*de*)-$
(19)
2 @(6*)fi(T’))x, de*
oi[xk(Ei)]
Jb
fi(T*)
Sin
r\kf7
-
r*)
dT*).
(20)
Substituting (20) into (9) and employing the result in (4) gives:
‘Uk, 7) =
kzzk-' (1.:’ (%I(~*) COShkr -t A/c-“&~*)
Sin&T
+At-’ $
OifXkb)]
l
f&f+)
Sin Ak(T-
T*)
dT*
>
- 2
Xk(E)
&*)fi(~)}x,
-I- $
gi(e)fi(T).
de* (21)
Solutionof time dependentboundaryvalue problemsby the boundary OPerarOr method
611
Noting the eigenfunction expansion: gi(E) = 2 4-l k=I
” gi(r*)Xk de* Xk(e), lr
(22)
(21) becomes: q(e, 7) = kz, &-’ (1:’ {q&*) COshk7
+Ak-’
$,
oi[xk(s)l
1:1’(T*)
+
hk-‘@O(e*)
Sin
hkT}Xk
de* (23)
sin hk(T - T*) dT*)xk(r)
for ~1Ce
+A,+-“&+*)Sin
$!(5 7) = & 4-l (j-y{*Ok*)
COShkT
+Ak-’ 2
I,’fi(T*)
oi[Xk(Ei)]
for e1 < E <
Sin
Ak(T
&.7)x,
- T*)dT*)
Xi(e)
de*
(24)
e2.
NUMERICALEXAMPLE The uniaxial stress is evaluated in a homogeneous bar for which the boundary conditions are: f
l(T)=~(0,7)=crosin(~)
fOiO’T~TcJ
=Ofor
T>Tq
iNf f2(7)=-(1,7)=0 ac Note that (25) implies that LY~ = a2 = +7r/2in (2). Thus the eigenfunctions in (9) satisfy free-free boundary conditions. From [7], the required eigenfunctions and eigenvalues for the homogeneous bar are cos(km) and kn, respectively. The initial conditions are taken as: Y&) =
*o(c)= 0.
(26)
From (23): w(C,7) = -2 kz, ck(r) T
I
O
(27)
where:
ck(T)
=
16fdT*) Sin ~T(T - T*) dr*.
(28)
The dimensionless stress is obtained from (24) as: g
(e, 7) = 2 kz, ck(r)
Sin(h),
0 <
E <
1
(29)
H. D. FISHER, M. M. CEPKAUSKAS and S. CHANDRA
612
With R, = T/T,, and p2 = kr, the C&(r) are given as follows for T I
_
sin& I
go COS =--2 (
,B 2(
-
PI--P2
PIT
Pt)r
_
sinU%
+ 82)~ I
pl+B2
cos
827
)
, p1 f
TO
sin 2&r - sin /3rr cos 2&7 + sin /?,r -7cos/3r1 > 2Pi
Cos8r7sin2/3rro -sin817cos2P170+sin8,7 26,
(30)
p2
PI =P2
-rocos817 >, fh=fh.
(31)
(33)
The behavior of (29) as E +O has been explained by Friedman[l3]. For an infinite number of terms in (29), this reference ~emons~ates that:
where ~(7) is determined from the expansion: ck(7)=
2 y.
(35)
m-l
Integration of (28) by parts results in c,(r) = fi(r)/n. Thus from (34):
For a compressive pulse with amplitude a0 = -1.00 and a time duration 70= 0.50, the digital computer program WAVE, available from the authors upon request, was employed to compute various displacement and stress profiles. Equation (27) was utilized together with the rigid body dispIa~ement, -lo’ (co(r - r*) sin (~~*~?o) dT* from Apjmd~ 2 Of [71, to ascertain the axial variation of bar displacement at T = 0.25, 0.50, 0.75 and 1.00, respectively. The results are shown in Fig. 1. Figures 2 and 3 present curves calculated from (29) for T = 0.25 and KS = 10,100and 1000, respectively. At this time the amplitude of the pulse at the origin is unity and the leading edge of the pulse has arrived at the location E = 0.25 [tco = (~T/CO~CO~ = LT, E = LTIL = ~1. Clearly inclusion of more terms in (29) improves the sedation of the actual stress crtndition by decreasing the size of &heregion where the pulse rises from a spurious value of zero to the required unit value of fi(0.25) and by eliminating the extraneous oscillations near the leading edge of the wave, Fig. 2 also demonstrates that the marked improvement in rise time of the plotted solution for KS = 1000 compared to the curve for K, = 100 is associated with negligible increase in overshooting since the first stress peaks are identity to plotting accuracy. Figure 4 is a plot of two stress profiles in the time regime after pulse application. Since reflection from the free end changes the sign of the incoming stress wave, the spatial variations of the stress at r = 1.5 and T = 2.0 are equal to the negative of the solutions obtained for T = 1.0 and 0.5, respectively. SUMMARY
Review of the above numerics exampie discloses that in order to be assured of the validity of the solution of interest, either in other second order systems or in higher order systems such
Soiution of
timedependent boundary v&e problems by the boundary operator method
613
as gCw~f% beams azxl p&s, *43 cond%Qas mw& z3@ msa. T&se s%z &at ti ~~~~~ ~~~~~~~~ d zttztnyseries to Zs ~~~~~~~~ txttts lx prwen (fhis was &tte for (24) in {143and
that comparable behavior to that of (36),must be shown at the boundaries of the region which are subjected to time-dependent excitation, In problems where these requirements are satisfied, the analyst gains three essantial benefits, in addition to the relative simplicity af the mathematics, from ern~l~~i~ the IX?M. First, the solutions for the de~e~deot variables ttf interest are not ~s~cted to a limited range of the temporal variable or a particular set of anon conditions. In the problem investigated above, the effect of I&-i&let (fixed), Neumann ffr~e) or Cauchy (end constrained by a light spri~~~l5]~,~undary perditions can be ascertained by setting ~1% ty2= 0, (~1 = (~2= n/2, aI, (~2P 0 f 42, respectively in (2). Second,
c
-1.0
b -0.9
0
0.001
o.am
o.w3
o.wd am0.m o.wd DIMENSIONLESS POSITION, (
0.m
0.039
Fig. 2, Stress solution near (0 c: r ~0.0101) forced boundary.
Fig. 3. slrr?lastoiution in vicinity (0.01CC 6 I 0.50) of forced boundary.
0.010
614
H. D. FISHW,hf. hf. CEPKAUSKAS and S. CHANDRA
0.2
0.3
0.4
0.5
DIMENSIONLESS
0.6 POSITION
0.7
0.8
0.9
1.0
E
Fig. 4. Stress solution after pulse application.
since auxiliary functions are not employed, the uniqueness of the solution is insured. Third, the necessity of determining and evaluating particular auxiliary functions corresponding to designated boundary conditions which is a requirement of [IO]has been eliminated. REFERENCES I. G. Cinelli, Dynamic vibrations and stresses in elastic cylinders and spheres. J. Appl. Mech. 33, 825-830(1966). 2. T. C. Lee, Note on wave propagation in a finite nonhomogeneous rod. J. Appl. hfech. 41,291-293 (1974). 3. H. L. Schreyer, Inverse solutions for one-dimensional seismic waves in elastic, inhomogeneous media. 1. Appl. Me&. 44,469-414 (1977). 4. S. C. Chou and R. Grief, Numerical solution of stress waves in layered media. AIAA J. 6(6), 1067-1074(196@. 5. R. Grief and S. C. Chou, The propagation of radially symmetric stress waves in anisotropic nonhomogeneous elastic media. 1. Appl. hi& 38, 51-57 (1971). 6. W. E. Alzheimer and hf. J. Forrestal, Elastic waves in thick-walled cylinders, 1. Engng Mech. Div., AXE %(EM2), 193-199(1970). 7. U. S. Lindhohn and K. D. Doshi, Wave propagation in an elastic nonhomogeneous bar of finite length. 1. Appl. Mech., 135-142(1%5). 8. H. D. Fisher, Wave propagation in a hollow conical frustum. J. Appl. Mech. 39, 1159-l161(1972). 9. Y. H. Pao and A. N. Ceranoglu, Determination of transient responses of a thick-walled spherical shell by the ray theory. 1. Appl. Mech. 45, 114-122(1978). IO. R. D. Mindlin and I. E. Goodman, Beam vibrations with time dependent boundary conditions, 1. Appl. Mech. 17, 377-380(1950). 11. L. Meirovitch, Analytical Methods in Vibrations,pp. 300-308.Collier-Macmillan Canada, Toronto, Ontario, (1%9). 12. H. I. Epstein, Qrations with timedependent internal conditions. J. Sound and Vibration,B(3), 297-303(1975). 13. B. Friedman, Principles and Techniques of Applied Mathematics, p. 271. Wiley, New York (1956). 14. W. Kaplan, Advanced-Calculus,p. 653. Addison-Wesley, Reading, Massachusetts (1959). 15. C. C. Lin and L. A. Segal, MathematicsAppliedto DeterministicProblems in the Natural Sciences, p. 372.Macmillan.New York (1974).