FINAL REPORT 

on 

FATIGUE BEHAVIOR OF REXHOSES & BELLOWS WE TO ROW-INDUCED VIBRATIONS 

Sponsor: 

NASA Kennedy Space Center 

Sponsor Contract No: NAG10-0017 

Prepared by 

Prateen V. Desai 

Pri nci pal Investigator 

and 

Lindsey Thornhi 1 1 

Graduate Research Ass? stant 

George W. Woodruff School of Mechanical Engi neeri ng 

Georgia Institute of Techno1 ogy 

January 1986 
ABSTRACT 

This report summarizes the analysis and results developed in a fresh 

approach to calculate flow induced vibration response of a flexible flow 

passage. The vibration results are further examined in the frequency domain 

to obtain dominant frequency information. A cumulative damage analysi s due to 

cyclic strains is performed to obtain number of cycles to failure for metallic 

be1 1 ows of particular speci fications under a variety of operational 

conditions. Sample plots of time and frequency domain responses are included 

in Appendix I. complete listing of a computer program is provided in this 

report as Appendix 11. The program successively executes each of the analyses 

needed to calculate the vibration response, the frequency response, the cyclic 

strains and the number of cycles to failure. The program prompts the user for 

necessary input information. Sample data from the program is provided in 

Appendix 111. The fatigue life results obtained by the computer model lie 

within an acceptable range of previously measured avail able data. 
I. INTRODUCTION 

Project Background and Overvi ew 

This report describes all the work performed by Georgia Institute of 

Techno1 ogy, Col lege of Engi neeri ng , School of Mechanical Engi neeri ng under 

Grant number NAG10-0017, "Fatigue Behavior of Flexhoses & Bellows Due to .Flow 

Induced Vibrations". This study was performed for the Kennedy Space Center of 

the National Aeronautics and Space Administration. 

The general objective of this work was to investigate the physical 

phenomena associated with flow-induced vibrations of a flexible line in order 

to develop a methodology to examine its fatigue li fe. In response to NASA 

needs, much emphasis was placed on developing a computational procedure for 

calculating the number of cycles to failure for a speci fic be1 lows based on 

the assumed vibration model. 

Flexible expansion joints are currently being used in a variety of 

applications. These include the nuclear reactor cooling systems, supply and 

return loading lines for submarines from the mother ship on the base, liquid 

cryogenic fuel lines for the space shuttle external tanks and so forth. There 

i s, in fact, a manufacturer's association for expansion joints which pub1 i shes 

its recommended practices for their use. 

Serious di fficul ties are encountered when practical solutions are not 

obtai ned to control the so-cal 1 ed " sel f-generated" sound or vibrations in fl ow 

systems which use such corrugated lines. For the special case of the 

internal ly corrugated flexible hoses, segements of which are utilized in the 

feed lines of Space Shuttle external tanks, a qualitative mechanism was 

identified in terms of the fluid flow related parameters governing the 

physical phenomena 1. Unshrouded hose vibration was attributed to an 

approximate matching between the longitudinal structural frequency and the 
frequency of vortex shedding at a critical flow velocity over the first few 

corrugations in the hose. This led to the computation of the modal 

frequencies of vibration by representing the hose with a series of lumped 

spri ng-mass elements and re1 ati ng the ampl itude of vibrations to the damping 

characteristics of the system. It was found necessary to supplement the 

postulated model by a large number of empirical correlations obtained from a 

series of experiments to relate the system geometry, the mean flow rate and 

the fl uid properties with the modal frequenci es and ampl i tudes of vibration, 

the damping characteri stics, the maximum a1 ternating stress and the expected 

fatigue li fe of the hose material [Z]. 

Experimental observations under simul ated fie1 d conditions have proved 

the predicted failure results of the two cited works to be very 

conservative. Furthermore, an overwhelming re1 i ance on empirical ly obtained 

functions to predict the fagitue life of the hose in these studies has made it 

cumbersome to di rectly uti 1 i ze the suggested fai 1 ure prediction procedure. 

Flex line flow-induced vibrations were a1 so examined in another study [3] with 

an eye to developing a methodology to assess its fatigue life. The cavity 

resonance model utilized in this work tends to be extra sensitive to minor 

parameter variations in predicting the time to failure by fatigue stresses. 

It appears that a definite need exists in developing a vibrational mode1 for 

the flexible line that incorporates the interactive dynamics between fluid 

flow past cavities and the flow-induced vibration of the convolutes of the 

be1 1 ows configuration. 

In the context of flow past cavities, there are three types of self 

generated oscillations. These are the fluid-dynamic, the fluid-resonant and 

the fluid elastic types of oscillations [4,5]. The fluid dynamic oscillations 

arise from some type of flow instability inherent to the flow configuration. 
The fl ui d-resonant osci 11 ati ons occur due to standi ng wave resonance effects 

which are particularly important for compressible fluid flow. The fluidelastic oscillations involve a coupling of the sol id boundary vibrations with 

the fluid flow perturbations causing those. It is this coupling which must be 

quantitatively incl uded in the vibration model which predicts the convol ute 

deformations. Subsequent cycl ic strain analysi s and the fatigue damage due to 

it are a1 so examined in thi s work. 

Statement of Research Problem 

This research concerns the development of a vibration model as well as 

fatigue damage analysis to estimate life cycles of a bellows or a flexhose 

with a single- or multiply wall construction. The vibration model emphasises 

the coupled interactive nature of the fluid flow and structural vibrations. 

The vibration response is further examined to assess the dominant frequencies 

causi ng fatigue damage. Fi nal ly, the research adopts exi sti ng theories of 

cumulative fatigue damage with the results of vibration analysis to assess the 

number of cycles to failure of the be1 lows. 

Flow Past Cavities 

One of the first modi fications to improve currently existing models of 

flow excitation in a bellows geometry concerns the force acting on each 

individual convol ute. A preponderence of previously pub1 i shed information 

concerns flow past rectangular cavities. Of particular interest is the 

pressure distribution along the walls of the cavity. 

Sketches of flow patterns in a rectangular cavity based on photographs 

[6] are shown in Figures l(a), (b) and (c). 

A plot of the measured pressure distribution [6] along the cavity wall 

for three different values of aspect ratio, h/b (0.5, 2.0, 1.5) reveals a 

trend shown in Figure 2. Similar plot [7] is shown in Figure 3. On the other 
Figure 1. Sketches of Flow Past Square Cavities Based on 

Photography [6]. 

Figure 2. Typical Measured Pressure Distributions Along Cavity 

Wall Showing Influence of Height-to-width Ratio 16'1. 
1 

Left Slde Bottom Right Side 

(B Experiment 

- Predlctlor 

. . 

Figure 3. Pressure Distribution Along Square. Cavity Walls ['I]. 

- CLVfTY 

I. 

Physicol Row model 

Typuol longitudinal pressure dirtribution in rheor la ye^ 

0) 

.- 

Figure 4. Pressure Distribution From Analytical Model [6 1. 
hand, an analytical model and its results of pressure variation along 

horizontal axis [61 are shown in Figure 4. It is seen that there is a 

pressure maximum at an area of flow reattachment on the downstream wall of the 

cavity. 

On the assumption that these same trends hold for flow past a convoluteshaped cavity, the expected flow patterns are as shown in Figure 5. The 

maximum wall presure is shown acting near the top of the downstream wall as 

demonstrated in the earlier results on cavity flow. The drop in pressure just 

previous to the sharp pressure rise shown in the plots might be explained by 

the vortices in the recirculating flow on the downstream wall. The vibrations 

of the bellows cause these recirculating vortices to move along the wall and 

be shed from the convolute tip periodically. As a vortex crosses the region 

of maximum pressure it causes a drop in this pressure. After the vortex has 

been shed the maximum pressure is reestabl i shed, thus produci ng a fluctuating 

force. Thi s fl uctuati ng pressure force, along with an acoustical resonance 

( radi a1 ) that may exi st, could produce re1 ati vely high pressure level s 

resulting in a bending of the convolute. 

The maximum pressure force pictured may be broken down into transverse 

and axial components. Previous research has described the force acting on the 

convolute as a longitudinal force. Actual ly, thi s longitudinal force is seen 

in the present research as only the horizontal component of the maximum 

pressure force. 

11. DEVELOPMENT OF THE MATHEMATICAL MODEL 

Assessment of fatigue failure of metallic bellows due to flow induced 

vibrations requires 

a. the development of a coupled fl uid-structure interactive vibration 
Figure 5. Impact Pressure Due to Flow Past a Convoluted Cavity. 
model to predict individual convol ute deformations, 

b. a frequency domain analysis of the vibrational response to identify 

domi nant frequenci es. causing the damage 

and 

c. a failure analysis based on available cummulative damage theories to 

calculate the number of cycles to failure. 

Flow Induced Vibration Model 

The interactive dynamics of bellows vibrations due to fluid flow over a 

cascade of convoluted cavities, which in turn modi fies the fluid flow and 

forces that excite the vibrations, is a coupled nonlinear phenomena. It is 

governed by fluid flow conditions, the bellows geometry and the material 

response behavior. 

The effective excitation force is internally generated and an analysis 

must be performed to express it in terms of identifiable physical 

quantities. The vibration of an individual convolute within a mu1 ticonvol uted 

bellows in response to the excitation may be examined via lumped springs and 

masses in terms of physical quantities. 

The Excitation Force 

A comprehensive analysis of flow induced vibration of a bellows must 

consider the motion of an individual convol ute in two coordinate directions. 

These di rections correspond to the longi tudi nal and the transverse motion of 

the convolute as illustrated in Figure 6(a). 

An expression for the excitation force acting on the convolute must 

account for the variation of the force due to convolute motion. When the 

bellows are absolutely motionless the flow geometry is that illustrated in 

Figure 6(b). This flow may be represented by the velocity triangle in Figure 

7b). 

In this figure, a, , estimated from data on flow past cylinders [8],is 
Figure 6a. Coordinate System to. 

Analyze Convolute 

Deformation. 

Cigure 6b. Configuration for Unde-- 

formed Be1 1 ows . 

Figure 6c. Convolute Deformation Along x. 

... 

Figure 6d. Convolute Deformation Along y. 
the angle between the velocity of the flow approaching the convolute tip, Uo, 

and the velocity of the core flow, Uc. The flow configuration illustrated in 

Figure 6(b) produces force coefficients on the convolute tip given by 

and 

where CD and CL are drag and 11 ft coefficients, respectively, of convolute tip 

in the core flow; CD and CL are estimated from flow past cylinders [8]. 

When the convolute tip is di splaced longitudinally, the configuration 

illustrated in Figure 6(c) results, The angle made by the flow approaching 

the convolute tip with the core velocity, ai (= a o + ba) is shown on the 

velocity triangle in Figure 7(b) representing thi s situation. In this 

figure, Ur,x,i is the velocity of the flow approaching the ith convolute tip 

relative to the convolute longitudinal velocity. The force coefficients for 

this situation become 

and 

uL r,x ,i (C sin ai + CL cos ai) . cy,i=~ D 

up 

These expressions show that longitudinal motion of the convolute produces a 

vari ation in both the longitudinal and transverse excitation forces. 

Therefore, longi tudi nal motion cannot exi st independent of transverse motion, 

When the convolute tip is displaced in the transverse direction, the 
Figure 7a. Velocity Triangle: Undef ormed Be1 1 ows . 

Figure 7b. Velocity Triangle: Deformation Along x. 

u, 

Figure 7c. Velocity Triangle: Deformation Along y. 

Figure 8. Velocity Triangle Combined x & y Deformation. 
configuration shown in Figure 6(d) represents the situation. The new 

angle, made by the approaching flow with the core velocity (@i = a. + A$) 

is shown on the velocity triangle in Figure 7(c). In this figure Uryyyi is 

the velocity of the flow approaching the ith convolute tip relative to the 

convolute transverse velocity. The force coefficients for thi s variation are 

and 

Similar to the case for longitudinal motion, transverse motion causes 

variation in both the longitudinal and transverse forces acting on the 

convolute tip. Therefore, transverse motion cannot exi st independent of 

longitudinal motion. 

I f the transverse and longitudinal di spl acements are a1 lowed to occur 

simultaneously, the velocity triangle becomes that shown in Figure 8. From 

this figure the angle (i and the relative velocity UrSi may, respectively, be 

determi ned as 

and 

The force coefficients are given by 
and 

The excitation forces may, therefore, be expressed as 

and 

where p is the fluid density and A is the convolute tip surface area. 

Upon substituion of the expression for into the expressions for the 

force coefficients one obtains 

and 

u r,i Uctana0 + (4.j - Yi-1 

Cy,i = 7 ( cLcos [tan-'( )I 

Uc uc - - ii ) 

Upon introducing the trigonometric identities 

and 
the force coefficients, with some rearrangi ng , become 

and 

Use of Equation (5) into Equations (9a) and (9b) yields the desired 

expressions for the force coefficients as C93 

and 
The excitation forces may now be written as 

and 

It may be noted that in addition to fluctuations in the excitation force due 

to convolute motion, there is a1 so a fluctuation exressed quantitatively due 

to vortex shedding that was suggested in a qualitative way previously by 

others. As the vortex moves across the region of highest pressure on the 

convolute tip, the pressure is reduced. After the vortex has been shed the 

hi gh pressure region i s reestabl i shed, thus producing an "on-off-on" type of 

fluctuation in the force which is already varying with convolute motion. 

Equations (lla) and (llb) are to be incorporated as excitation forces in 

the vibration model that represents the dynamic force balance of the be1 lows. 

The Lumped Parameter Vibration Model 

The development of an expression for the excitation force acting on an 

individual convolute suggests the need for a mechanical model which a1 1 ows 

both longitudinal and transverse motion of the convolute. In order to develop 

such a model the bellows is divided into discrete mass elements such as the 
one shown in Figure 9. By considering longitudinal and transverse motion of 

this mass element and the elements adjacent to it, a two dimensional vibration 

model is developed. This is illustrated in Figure 10, where FD,x and FD,~ are 

damping forces which represent such effects as sliding between plys of the 

bellows and internal frictional damping. Viscous damping is not shown in the 

illustration, since the damping effect produced by the fluid surrounding the 

convolute is a1 ready included in the expression for the excitation. kx and ky 

are material spring constants, rn, the material mass of one element and mf the 

fluid mass of one element. These components of the lumped model are further 

defined el sewhere in this report. 

From the vibration model illustrated in Figure 10 a system of equations 

describing the motion of the convolutions must be developed. This system must 

include a1 1 of the lumped parameters and effects described above. The system, 

when solved, must a1 so provide displacement as a function of time for each of 

the convolute tips-. 

The development of the equations governing the system behavior results in 

two times the number of convolutions coupled, nonlinear, ordinary, 

di fferenti a1 equations. These equations are derived by performing a force 

balance on the discrete mass element shown in Figure ll(a). The force balance 

gives 

which upon rearranging becomes 
Figure 9: Convolute Representing One Discrete Mass. 
Figure 10. Vibration Model. 

- 19 - 
Figure Ila. Force Balance Along x-direction. 

Figure Ilb. Force Balance Along y-direction. 
The force balance in the transverse direction also may be performed on the 

element shown in Figure ll(b). This gives 

Once agai n , rearrangement gi ves 

Upon substitution of the previously determined expressions for Fx,i and Fy,i 

the system governing equations become. 

and 

.. ('"m + mf)yj - FD,y sgn(ii) + Zkyyj - ky~i+l - 'y~i-l 

- PA {CU, tanao + (ii - il~2 + [uC - (ii-l 2 112 

- 7 - ii)l 1 

{CLCUc - - if 11+ CDCUCtana, + (4i - yiSl )I3 , 

where each value of i designates a speci fic convolute. This produces the 

previously mentioned system of coupled, nonlinear, ordinary, differential 
equations. 

The solution of this system is sought under the assumption that the 

bellows is rigidly attached at each end and that the entire system is 

i ni ti a1 ly quiescent. 

Components of the Lumped Parameter Model 

The Fluid Added Mass 

As a convolute is displaced it also displaces a certain amount of the 

fluid flowing through the bellows. This implies that the mass of this fluid 

must be included in the governing equations as an inertial effect; hence, the 

.. .. mf xi and mf yi terms. Each convolution has an instantaneous configuration 

ranging between the two extremes shown in Figure 12(a). This suggests that 

the fluid added mass varies with instantaneous longitudinal position much more 

than transverse position. The expression for the fluid added mass will 

therefore be a function of the longitudinal position of adjacent convolutes. 

With the geometry defined in Figure 12(b) an expression may be developed for 

the fluid mass. For example, the instantaneous cavity width and length of the 

straight section are given by 

and 

With this information the areas A1, A2, and A3 may be determined and swept 

around the bellows to yield the volume of the fluid added mass. Multiplying 

this volume by the fluid density, p , yields the expression for the fluid 

added mass, mf, as 
Figure 12a. Instantaneous Fluid Mass Configurations, 

Figure 12b, Geometry to Compute Fluid Mass. 
Here, Dm is the mean bellows diameter. Since this expression is multiplied by 

.. .. xi and y in the governing equations it introduces additional nonlinear i 

effects into the system of governing equations. 

The Materi a1 Mass 

Each discrete mass in the vibration model is composed of the material 

mass of one convolution as well as the fluid added mass. Having determined an 

expression for the fluid mass added to one convolute in m6tion, an evaluation 

of the materi a1 mass which composes one convol ute may a1 so be obtained. 

A cross-section of one convolute is illustrated in Figure 13. The area 

of this cross-section is given by 

AREA = [(2r + 4)a + 2hlt (19 

The volume of material may then be determined by sweeping this area around the 

be1 1 ows mean diameter to yi el d 

VOLUME = nDm[(2n - 4)a + 2h3t . (20) 

The material mass of one convolute, %, is determined by multiplying this 

volume by the material density, pm . This gives 
Figure 13. 

Figure 14a. Configuration to Calculate 

Spring Constant, kx. 

Convolute Geometry 

to Calculate Material 

Mass. 

Figure 14b. Configuration to 

Calculate Spring 

Constant, k Y' 

Figure 15b. Geometry to Calculate 

Figure 15a. Geometry to Calculate - 25 -. Convolute Tip Strain 

Convolute Tip Strain Along y-direction. 

Along x-direction. 
The Longi tudi nal Spri ng Constant, k, 

The longitudinal spring force is the force tending to restore the 

convolute to its equilibrium position when it has been deformed by a 

longitudinal displacement. This force is a result of the elasticity of the 

be1 lows materi a1 . In the vibration model the spring force is represented as a 

longitudinal spri ng constant, kx, mu1 tip1 i ed by the di spl acement of the 

convol ute re1 ati ve to adjacent convol utes. 

An expression for the spring constant can be determined by using the 

strain energy method along with the theorem of Castiglipno. From the geometry 

illustrated in Figure 14(a) and with a force applied in the longitudinal 

direction the strain energy of this configuration is given by 

Here, U is the strain energy for the configuration illustrated in Figure 

14(a), E is the bellows material modulus of elasticity, and G is the shear 

modulus. A1, 11, are the cross-sectional area and moment of inertia for the 

section based on the bellows inside diameter. A2, 12 and A3, I3 are based on 

the be1 lows mean diameter and outside diameter, respectively. 

The theorem of Castigliano states that the derivative of the strain 

energy with respect to an applied force is the displacement of the point of 

application of the force in the direction of the force. Therefore, the 

longitudinal displacement of the convolute tip illustrated in Figure 14(a) may 

be determined as 
Since the longitudinal spring constant is the fdrce required to produce a 

given di splacement, the spring constant may be determined as 

The Transverse Sprl ng Constant, ky 

The transverse spring constant may be determined in a manner anlagous to 

the method used for determining the longitudinal spring constant. The strain 

energy for the configuration illustrated in Figure 14(b) is found to be 

By the theorem of Catigl iano the transverse di splacement is given as 
The transverse spring constant may now be determined as 

111, SOLUTIONS IN THE TIME AND FREQUENCY DOMAINS 

Solution of the Governing Equations In the Time Domain 

The bellows vibrational response to flow excitation is obtained by 

solving the system of governing equations for the displacements. Since an 

analytical solution to this system is not known a numerical solution to be 

executed by computer must be chosen. In problems where the range of 

integration i s considerable, the numerical method empl oyed must be stab1 e 

and/or relatively stable to obtain an accurate solution. A1 so, in a vibration 

problem where the motion is to be studied over a substantial number of 

oscillations, a method having a ma1 1 per-step truncation error should be used 

in order to minimize the cumulative error. 

Upon consideration of these requirements, a "predictor - modifier - 

corrector" method called Hamming's Method was selected as the solution 

technique [lo]. Hamming's Method is both stable and relatively stable and 

results in small per-step truncation errors. 

The iteration algorithm for a di fferential equation of the form 
may be stated as the predicted value of y at t + 1: 

P(yt+1) = Yt-3 

4h 

+ J C2YIt - ~'~-1 + 2~'~-~ I 

The modi fi ed val ue 

11 2 M(Y~+~) = P(Y~+~) - [P(yt) - C(Yt)] 

The corrected value 

1 C(Y~+~) ' 3 {9~t - Yt-2 + 3h[m(~'~+~) + 2~'~ - ( 28d 

The final value 

9 

F(yt+1 = C(Yt+1) + CP(Yt+,) - C(Y~+~)~. 

Here, h is the time step used in the iteration procedure. 

As can be seen from these equations Hamming's Method a1 so has the 

advantage of being relatively efficient in terms of computing time, since 

there is only one function evaluation per time step. A di sadvantage, however, 

is that Hamming' s Method is not sel f-starti ng. It requires solution values at 

the first three time steps in order to begin the iteration process. 

To provide the necessary solution for the first three time steps, a 

fourth-order Runge-Kutta procedure i s employed. While being very inefficient 

in terms of computing time, the Runge-Kutta procedure does have the advantage 

of being self-starting. 
By starting the iteration with an initially quiescent system, the RungeKutta procedure determines the longitudinal and transverse di spl acements of 

each convolute tip for the first three time steps, This information is then 

passed on to Hamming's Method which iterates to 4096 time steps, each step 

being of 0.01 mil 1 i second duration. The di splacement information for the last 

2048 time steps is stored as a "data sample" for frequency and fatigue 

calculations. By choosing the last 2048 time steps the "data sample" is 

considered far enough from time zero where steady-state vibration is assumed. 

With the 0.02048 second sample of steady-state vibration di splacement 

information, a frequency response analysis and strain analysis may be 

performed. Plots of displacement versus time from this solution are shown in 

Appendix 1. 

Frequency Response 

Examination of the displacement versus time plots reveals that both 

longitudinal and transverse vibrating of an individual convolute are made up 

of several frequency components. Each of these frequency components 

contribute to the damage done to the convolute during a given time period. In 

order to estimate the extent of this damage on each convolute the dominating 

frequencies of vibration in both the longitudinal and transverse directions 

must be determined. 

With the time-dependent di splacement obtained from the solution of the 

governi ng equations, a Fourier Analysi s may be performed. Thi s analysi s 

yields the frequency response of each convolute in each coordinate 

direction. The Fourier Analysis is performed by means of the Fast Fourier 

Trans from a1 gori thm which i s developed here. 

By representing the convolute response as a complex Fourier series one 

obtains 
and 

where X(U) = 1- x(t)e-i't dt 

w is a discrete angular frequency and t is time. Recall from the solution of n 

the governing equations that the sampling period, T, and the number of data 

samples, N, are respectively, 

T = 0.02048 second 

N = 2048 samples 

Therefore the time step size is 

The fundamental frequency, U, , is given by 
and the maximum frequency which may abe determined, omax , is 

W 

- 2s - 2s -7- , Noo . max ~t Im 

Define a discrete set of N frequencies at which the convolute response 

wi 1 1 be dete rmi ned. 

Since x(t) and y(t) vanish for time less than zero and since they have not 

been determined for time greater than T, the transform may be approximated as 

and 

If these integrals are approximated by summations, the expressions become 

and 

Upon evl uating these transform expressions at di screte frequencies, W, , one 

obtains 
Let WN = exp[-i $ 1 to get 

and 

Each power of WN increases the polar angle by 2r/N . After N powers the 

pattern repeats so it is not necessary to compute all powers. Therefore, 

define Uk such that 

and 

where MOD(k,N) is the modulus function. The combination of Uk with the 

previous transform expression yields the equations for the transformed 

response at di screte frequenci es as 

and 
Note that since (Wn) m(-n) is the complex conjugate of (w~)'~, it is sufficient 

to evaluate X(wn) and Y(on) for n , 0 . 

The results of the Fast Fourier Transform are the frequency response 

plots shown in Appendix 1. With this information the dominant frequencies of 

both longitudinal and transverse vibration of each convolute may be 

determined. These dominant frequencies will later be used to assess the 

cumulative damage effects on each convolute, 

IV. CYCLIC STRAIN, DAMAGE AND FATIGUE LIFE 

Estimation of Cyclic Strains 

The numerical model to calculate fatigue life of bellows is based on 

obtaining strain amp1 i tude for the cycl ic deformation of each convol ute. The 

strain amplitude has to be evaluated at the convolute tip since this is the 

location at which the bellows is known to fail. A maximum and minimum value 

for the tip strain resulting from both longitudinal and transverse vibration 

needs to be determined, Since the strain induced at a point by di fferent 

mechani sms is additive (same fibres damaged), the maximum tip strain resulting 

from 1 ongi tudi nal vibration and the maximum tip strain resulting from 

transverse vibration may be added to get a total maximum convolute tip 

strain. Likewise, the corresponding minimum strains may be added to get a 

total minimum convolute tip strain, The strain amplitude may then be 

determined by the equation: 

E max AE = - 'min 

2 

where E is the tip strain and AE the strain amplitude. 
In order to determine an expression for the convolute tip strains 

resul ti ng from both 1 ongi tudi nal and transverse di spl acements the strai n 

energy method along with the theorem of Castigliano is used. For longitudinal 

displacement of the convolute tip it is assumed that the strain energy 

resulting from the displacement is concentrated in the section of the 

convol ute ill ustrated in Figure 15(a). The strain energy for this 

configuration is given by 

where, as before, E is the bellows material modulus of elasticity and G is the 

shear modulus. A1 and I1 are the cross-sectional area and the moment of 

inertia based on the be1 lows inside diameter. A2, I2 and A3, I3 are based on 

the be1 1 ows mean diameter and out side diameter, respectively. 

By applying the theorem of Castigl i ano the longitudinal di splacement of 

the tip, xi, may be determined as 

Since = EE x ,i , the displacement may a1 so be expressed 
and the convolute tip strain may be written as 

Since the transverse displacements are smaller than the longitudinal 

displacements the strain energy resulting from transverse di spl acements is 

considered to be concentrated in a smaller section of the convolute, Figure 

15(b) is an illustration of this section, Once again the strain energy for 

this configuration is given by 

and application of Castigliano' s theorem gives 

Upon substituting a = E . E and solving the resulting expression for E Y,i Y,l Y,i ' 

the equation for the convolute tip strain resulting from transverse vibration 

is obtained as 
As described earlier, the maximum values of E,,~ and E which occur during Y ,i 

the sampling period are added to get a total maximum tip strain, oi . The 

minimum values of E~,~ and E are a1 so added to yield a total minimum tip Y,? 

strain, E,~ . The convolute tip strain amplitude, Ari , may then be 

determined by 

Estimate of Cycles to Fai 1 ure 

A modern approach to characterize the fatigue behavior of material s is to 

focus on the cyclic strain li fe [12]. The total strain, E , is considered as 

having an elastic re , and a plastic, E components. Expressed as strain P' 

amp1 i tudes, thi s imp1 i es 

The elastic strain-li fe can be expressed as 

I 

where ua and of , respectively, represent the true stress amplitude and the 

fatigue strength coefficient, E is the modulus of elasticity, b is the fatigue 

strength exponent (typically, b varies between -0.05 and -0.12) and Nf is the 

number of cycles to failure. The plastic strain-li fe is related by a power 

function as 
I 

where cf is the fatigue ductility coefficient and d is the fatigue ductility 

exponent (typically, d varies between -0.5 and -0.7). The total strain 

amplitude may thus be expressed as 

This equation is cal led the strain-li fe relation and forms the basis for the 

strai n-1 i fe approach to predicting fatigue behavior of such materi a1 as 

wrought metal s. For thin-wal led she1 1 s, such as the be1 lows, Nf is typically 

taken to be the number of cycles to initiate a crack. 

With the strain amplitude for each convolute in hand it is possible to 

compute the number of cycles to failure for each convolute. For the types of 

steels from which be1 lows are commonly made, the equation constants are given 

of' = 150 kpsi 

and d - -0.7 

The strai n-1 i fe relation obviously requires an iterative solution. To 

accompl i sh thi s, simp1 e Newton-Raphson Method has been employed. Since the 

slope of the strain-li fe relation is always negative, there should be no 
problems with convergence so long as Nf is initially relatively mall in the 

iteration process. 

Cumulative Damage and Fatigue Life of Bellows 

By estimating the cumrnulative damage done to each convolute during the 

sample period an estimation can be made as to which convolute will fail 

first. The assumption made here is that the convolute undergoes cyclic strain 

of amplitude Asi for the time of the sampling period at each of the dominant 

frequencies determined by the Fourier Analysi s. For example: 

o convolute #3 vibrates with strain amplitude Ac3 at a frequency fl 

for a time equal to the sampling period. 

o convolute #3 then vibrates with strain amplitude As3 at frequency f2 

for a time equal to the sampling period 

o etc. 

From this approach the cumulative damage done to each convolute during 

one sampling period can be determined. This is done by use of the PalmgrenMiner linear-cumulative-damage rule [12] which may be stated as 

2N1 (Reversals applied at ual) 

eversal s to failure at ual) ' d=~=(~ f '1 

where d = damage. Failure will occur when 

In the present model the damage for each convolute during the sampling 

period may be determined by 
Here T is the sampling period, Nfi is the cycles to failure for the ith 

convolute, fi ,. .. fi ,, are the dominant frequencies of vibration, and di i s the 

damage for the fth convolute. Since steady state vibration is assumed to 

exist, the largest value of di corresponds to the convolute which is likely to 

fail first. 

The number of cycles to fai 1 ure for thi s convolute i s a1 so the cycles to 

failure for the entire bellows. A failure time wf ndow may be determined by 

dividing the bellows cycles to failure by the maximum and minimum dominant 

frequencies of the convolute which is first to fail. This defines a minimum 

time to failure as 

and a maximum time to failure as 

This tmi, and tmax then establish a failure time window. 
V. CLOSING REMARKS 

In testing the bellows computer program it became apparent that the best 

results are obtained for bellows with up to five convolutions. For bellows 

with greater than five convolutions the cycles to failure begins to decrease 

considerably. Therefore, it is recommended that program BELLOWS be used for 

bellows with up to five convolutions. 

The limitation of the program to give resonable vibration amplitudes for 

bellows with less than five convolutions is attributed to cumulative 

displacement effects in the springs connecting the discrete masses of the 

mechanical model. For more than five convolutions the displacements in each 

consecutive spring adds to eventual ly produce unreasonabl e convol ute tip 

di spl acement. Rather than introducing arti fi ci a1 damping into the governing 

equations to control these amplitudes, it was decided to limit the model to a 

five convolution model. Additional effects due to stress wave propagations in 

bellows with large number of convolutions need to be examined further in 

future research. 

For bellows with differing number of convolutions the interaction between 

the convolutes results in the strain energy being concentrated in a smaller or 

larger volume of the convolute. Fewer convolutions result in the strain 

energy being concentrated in a smaller volume near the convolute tip. To 

incorporate this effect into the model a factor multiplying the equations for 

strain c a l c u l a t i o n has been introduced. The strain multiplying factor varies 

linearly with the number of convolutions and its determination is internal to 

the computer program. 

A thorough literature search into the effects of sliding damping between 

plys and internal frictional damping did not yield sufficient information to 

explicitly incorporate these effects into the model. Just the same such 
effects do exist and were included in the model via damping forces. A 

comparison of calculated fatigue life cycles with measured values [2] indicate 

very encouraging agreements. The computer ca1 cul ations were made for 

equilvalent test conditions and bellows specifications to those used in the 

experiments, except that the number of convolutions used in the model was less 

than or equal to five. The agreement between the two indicated in Tables 1 

and 2 is within an order of magnitude. It is concluded therefore that within 

the scope of its application, the model presented in this study is quite 

sati sfactory . 

Table 1 Computer Results 

Be1 lows # Convol utes Flow Cycles Min. Time Max. Time 

Number in Model Velocity ( ft/s) to Fail (lo6) to Fail to Fail 

-( set) ( set 

- 

5028 5 3 5 2.9 979. 11,747. 

5028 5 50 210. 77,851. 856,362. 

5028 5 6 5 6.7 2482. 27,301. 

5013 3 3 5 0.67 651. 4556. 

5013 3 50 0.056 63.8 383. 

5005 3 35 0.33 285. 1708. 

5005 3 50 0.025 24.5 171.2 

Table 2. MSFC Tests 121 

Be1 1 ows Flow Cycles t g Time to 

Number Vel oci ty ( ft/ s)' Fail (10 ) Fail (sec) 

5028 65.8 ------ 787. 

5028 66.6 .I6 150. 

5013 5 0 1.3 2007 - 2027 

5013 45 -47 - -48 749 - 760 

5005 33.2 -78 - .79 1841 - 1863 

5005 65.2 -ma-- 75 - 100 

L 
VI . BIBLIOGRAPHY 

1. Johnson, J. E., Doffenbaugh, D. M., Astleford, W. J. and Gerlach, C. R., 

"Be1 lows Flow- Induced Vibrations" , Final Report, Contract No. NAS8-31994, 

Southwest Research Institute, October 1979. 

2. Gerlach, C. R., Tygiel ski, P. J., Smyly, H. M., "Bellows Flow-Induced 

Vibrations", Report number NASATM 82556, Marshall Space Flight Center, 

1983. 

3. Desai, P. V., Wu, X. F., "Flow-Induced Vibrations in Flexible Components 

of Space Shuttle Feed Systems", Project Number E25-605, 1984. 

Naudascher, E., "From Flow Stabi 1 ity to Flow-Induced Excitation ," Proc. 

ASCE., J. Hydraulics Div., Vol . 93, No. HY4, Proc. Paper 5336, July 1967. 

Rockwell , D., "Prediction of Osci 1 lation Frequencies for Unstable Flow 

Past Cavities", ASME J. Fluids Engineering, Vol . 99, 1977, pp. 294-300. 

Haugen, R. L., Dhanak, A. M., "Momentum Transfer in Turbulent Separated 

Flow Past a Rectangular Cavity", Journal of Applied Mechanics, Sept. 

1966. 

Yang, W. J., "Forced Convection Heat Transfer in Interrupted Compact 

Surfaces", ASME JSME Thermal Engineering Joint Conference, Vol . 3, pp. 

105-109, 1983. 

Eastop, T. D., Turner, J. R., "Air Flow Around Three Cyli nders at Various 

pi tch-to-di ameter Ratios for Both a Longitudinal and a Transverse 

Arrangement", Trans. Inst. of Chem. E., Vol. 60, 1982. 

Desai, P. V., Thornhill, L., "Fatigue Behavior of Flexhoses and Bellows 

Due to Flow-Induced Vibrations," Interim Report for NASA/KSC, NAG10=0017, 

July 1985. 

Mei rovi tch, Leonard, "Elements of Vibration Analysis" , McGraw-Hill Co., 

New York, 1975. 

James, M. L., Smith, G. M. and Wol ford, 3. C. Applied Numerical Methods 

for Digital Computation, Harper and Row, 1977. 

Socie, D. F., Mitchell, M. R. and Caul field, E. M., "Fundamentals of 

Modern Fatigue Analysi s" , Fracture Control Program Report #26, Col 1 ege of 

Engineering, University of I1 1 i noi s, January 1978. 

Blevins, Robert D., "Flow-Induced Vibrations" , Van Nostrand Reinhold Co., 

New York, 1977. 
APPENDIX I 

SAMPLE PLOTS FOR TINE AND FREQUENCY DOMAIN RESPONSES 

In the time response plots, the core velocity, Uc, is taken as 45 

ftlsec. The horizontal axis indicates time in seconds and the vertical axis 

indicates di spl acement in inch. x and y, respectively, indicate longitudinal 

and transverse vibrations displacements. Subscripts on x and y indicate the 

data for the particular convolution number. 

In the frequency response plots, the vertical axis indicates the Fouriertransformed di spl acements and the horizontal axi s incidates frequency in 

radi ans/second. The core velocity, Uc, for a1 l plots is 45 ftlsec. 















RAD. FREQ. 




APPENDIX I1 

LISTING OF PROGRAM BELLOWS 
PROGRAM BELLOWS (INPUT, OUTPUT, TAPE5=INPUT, TAPE^-OUTPUT) 

C 

C ..................................................... 

C * * 

C * BELLOWS FATIGUE PROGRAn * 

C * 1/28/86 JC 

C * * 

C * * 

C * * 

C * INVESTIGATOR: DR. P. V. DESAI * 

C * ASSISTANT: L. D. THORNHILL * 

C * * 

C ..................................................... 

C 

C THIS PROGRAM COMPUTES AN ESTIMATE OF THE FATIGUE 

C LIFE OF BELLOWS UNDER A VARIETY OF OPERATIONAL 

C CONDITIONS. 

C 

CALL ENTER 

CALL RUNGE 

CALL Iwm 

CALL STRAIN 

CALL FFT 

CALL CYCLE 

CALL DAMAGE 

CALL RESULT 

C 

STOP 

END 

C 

C 

C 

SUBROUTINE ENTER 

C 

C ENTER GEOMETRIC, MATERIAL AND FLOW PARAMETERS AND 

C CALCULATE EQUATION CONSTANTS. 

C 

COMMON/BLOK~/ PI,DELO,PIDEN,ALO,TE,EM,~C,~CT,N 

COMMON/BLOK~/ XK,YK,FDX,FDY,DENSA 

COMMON/BLOKS/ Al,E,BOTX2,BOTY2 

DIMENS ION FDX (9) , FDY (9) 

C 

PRINT*, 'ENTER CONVOLUTE PITCH IN INCHES.' 

READ*, PITCH 

PRINT*, 'ENTER CONVOLUTE HEIGHT IN INCHES.' 

READ*, HI ' 

PRINT*, 'ENTER PLY THICKNESS IN INCHES.' 

READ*, TPLY 

PRINT*, 'ENTER NUMBER OF PLYS.' 

READ*, NPLY 

PRINT*, 'ENTER INSIDE BELLOWS DIA. IN INCHES.' 

READ*, DIN 

PRINT*, 'ENTER NUMBER OF CONVOLUTES.' 

READ*, N 

PRINT*, ' ENTER FLUID DENSITY IN LBM/ (IN**3) . ' 

READ*, DENF 

PRINT*, 'ENTER FLOW VELOCITY IN FT/S.' 

READ*, UFEET 

PRINT*, 'ENTER THE ELASTIC MODULUS IN PSI.' 

READ*, ~l 

- 66 - 
PRINT*, ' ENTER THE MATERIAL DENSITY IN LBM/ (1~**3) . ' 

READ*, DENM 

C 

TE=TPLY*FLOAT (NPLY) 

UC=UFEET*12. 

E=E1*12. 

+~/(2.*(1.+1./3.)) 

C 

PRINT*, 'CONVOLUTE PITCH IS ',PITCH,' IN.' 

PRINT*, 'CONVOLUTE HEIGHT IS ',HI,' IN.' 

PRINT*, 'PLY THICKNESS IS ',TPLY,' IN.' 

PRINT*, 'NUMBER OF PLYS IS ',NPLY 

PRINT*, 'INSIDE BELLOWS DIA. IS ',DIM,' IN.' 

PRINT*, 'NUMBER OF CONVOLUTES IS ',N 

PRINT*, ' FLUID DENSITY IS , DENF, ' LBM/ IN**^) ' 

PRINT*, 'FLOW VELOCITY IS ',UFEET,' FT/S' 

PRINT*, 'ELASTIC MODULUS IS ',EI,' PSI' 

PRINT*, 'MATERIAL DENSITY IS ' ,DEW, ' LBM/ (1~**3) ' 

C 

DO 3 Jm1,N 

FDX (J) =8. 

FDY ( J) -2. 

3 CONTINUE 

C 

C CALCULATE EQUATION CONSTANTS 

C 

IF (uC .GT. 0. .AND. UC .LE. 180.) ALPHA-1.3 

IF (uC .GT. 180. .AND. UC .LE. 300.) ALPHA=1.2 

IF (UC .GT. 300. .AND. UC .LE. 420.) ALPHAII1.0 

IF (UC .GT. 420. .AND. UC .LE. 540.) ALPHA-0.6 

IF (UC .GT. 540. .AND. UC .LE. 660.) ALPHAo0.4 

IF (UC . GT. 660. .AND. UC . LE . 780. ) ALPHAmO .1 

IF (UC .GT. 780.) ALPHA=0.05 

R= PITCH/4. 

PI= ATAN( AN(^.) 

DOUT=DIN+~.*HI+~,*TE 

DM= (DIN+DOUT) /2. 

SIGMA- (2.*R)+TE 

DELO= (2. *R) -TE 

ALO= HI- (2. *R) 

PIDEN- PI*DM*DENF 

EM= PI*DM*DENM* (R* ( (2. *PI) -4. ) 4-2. *HI) *TE 

C 

AMOMl= (PI*DIN* (TE**3) /12. 

AMOM2= (PI-M* (TE**~. ) ) /12. 

AMOM3= (P I*DOUT* (TE**~. ) ) / 12. 

C 

AI=PI*DIN*TE 

A~=PI*DM*TE 

A~=PI*D~UT*TE 

B=HI-2. *R 

C 

BOTTOM= (~1/32. ) * ( (R**3 .) / (E*AMOM~) + R/ (E*A~) 

$ + R/ (G*A~) + R/ (E*A3) + R/ (G"A3)) + 0.125" (( (~""2,) / 

$ (E*AMOM~))*B + B/(E*A~)) + PI-16.)/32.)*((~**3.)/ 

$ (E"AMOM3)) 

BOTY2= (~1/32. ) * ( (R**3. ) / (E*AMOM~) +R/ (E*A~) +R/ (G*A~) ) 

C 

YK= l./BOTTOM 

C 

- 67 - 
RETURN 

END 

SUBROUTINE RUNGE 

CALCULATE STARTING VALUES FOR HAMMINGS METHOD. 

COMMON/BLOK~/ PI,DELO,PIDEN,ALO,TE,EM,UC,UCT,N . COMMON/BLOK~/ X,Y 

COMMON/BLOK~~/ XP,XPl,XP2,XP3,XP4,YP,YPl,YP2,YP3,YP4, 

$ XPP,XPPl,XPP2,XPP3,XPP4,YPP,YPPl,YPP2,YPP3,YPP4, 

$ CD, CL 

DIMENSION ~(0:6,0:2048) ,~(0:6,0:2048) ,xP(O:~) ,~~1(0:6), 

$ xP2(0:6) ,xP3(0:6) ,xp4(0:6) ,YP(O:~) ,~~1(0:6), 

$ YP~ (0:6) ,~~3(0:6) ,~~4(0:6) ,xPP(O:~) ,XPP1(0:6), 

$ XPP~ (0:6) ,XPP~ (0:6) ,xpp4(0:6) , YPP (0:6) , YPP~ (0~6) , $ ~PP2(0:6) ,yPP3(0:6) ,YPP4(0:6) ,cD(~) ,cL(~) 

COMMON/BLOK4/ TX,TXP,TY,TYP,K,Q,H,L,J,I 

REAL K 

DIMENSION TX (3) , TXP (3) , TY (3) , TYP (3) , K (4) , Q (4) 

TIME STEP, H, AND CONSTANTS. 

INITIAL CONDITIONS 
YP(J)=O.O 

YP1 (J)=O.O 

YP2 (J)=O.O 

YP3 (J)=O. 0 

YP4 (J) -0.0 

XPP1 (J)=O.O 

XPP2 (J) m0.0 

XPP3 (J) =O . 0 

XPP4 (J) -0.0 

YPPl (J)=O.O 

YPP2 ( J) -0.0 

YPP3 (J) 10.0 

YPP4 (J) =o. 0 

10 CONTINUE 

C. 

C BOUNDARY CONDITIONS 

C 

DO 30, JsO,N+l,N+l 

DO 20, I=O, 2048 

X(J,I)=O.O 

Y (J,I)=o.o 

20 CONTINUE 

30 CONTINUE 

C 

C RUNGE KUTTA STARTER. 

C 

DO 50, I=0,3 

c 

DO 31, J=l,N 

XP4 ( J) =XP3 (J) 

XP3 (J) =XP2 (J) 

XP2(J)=XP1 (J) 

XP1 (J)=XP(J) 

YP4 (J) =YP3 (J) 

YP3 (J) =YP2 (J) 

YP2 (J) =YP1 (J) 

YP1 (J)=YP (J) 

XPP4 (J) =XPP3 (J) 

XPP3 ( J) =XPP2 ( J) 

XPP2 (J) =XPP 1 (J) 

YPP4 (J) =YPP3 (J) 

YPP3 (J) =YPP2 (J) 

YPP2 (J) =YPPl (J) 

31 CONTINUE 

C 

DO 40, J=l,N 

C 

DO 32, M=1,3 

ID J+(2-M) 

TX (M) =X (IND, I) 

TXP (M) =XP 1 (IND) 

TY (M) =Y (IND, I) 

TYP (M) =YP1 (IND) 

32 CONTINUE 

C 

L=1 

CALL KQ 

C 

DO 34, M=1,3 

IW- J+ (2-M) 
TX (M) =X (IND, I) +A3*XP1 (IND) 

TXP(M)=XPl (IND)+ (~(1) /2.) 

TY (M) =Y (IND, I) +A3*YP1 (IND) 

TYP (M)=YP~ (IND)+ (~(1) /2.) 

CONTINUE 

L=2 

CALL KQ 

DO 36, M=1,3 

INP J+ (2-M) 

TX (M) =X(IND, I) +A~*xP~ (IND) +~4*K(1) 

TXP (M) =XP1 (IND) + (K (2) /2. ) 

TY (M) =Y (IND, I) +A3*YP1 (IND) +A4*Q (I) 

TYP (M) =YP1 (IND) + (Q (2) /2 .) 

CONTINUE 

L=3 

CALL KQ 

DO 38, M=1,3 

IND=J+ (2-M) 

TX (M) =x (IND, I) +H*XP 1 (IND) +U*K (2) 

TXP (M) =XP 1 (IND) +K (3) 

TY (M) =Y (IND, I) +H*YPl (IND) +A~*Q (2) 

TYP (M) =YPl (IND) +Q (3) 

CONTINUE 

L-4 

CALL KQ 

CONTINUE 

CONTINUE 

RETURN 

END 

SUBROUTINE KQ 

COMMON/BLOKI/ PI,DELO,PIDEN,ALO,TE,EM,~C,UCT,N 

COMMON/BLOK~/ XK,YK,FDX,FDY,DENSA 

DIMENS ION FDX (9) , FDY (9) 

COMMON/BLOK3/ X,Y 

COMMON/BLOK~~/ XP,XP~,XP~,XP~,XP~,YP,YP~,YP~,YP~,YP~, 

$ XPP,XPP1,XPP2,XPP3,XPP4,YPP,YPP1,YPP2,YPP3,YPP4, 

$ CD, CL 

DIMENSION ~(O:6,0:2048) ,~(0:6,0:2048) ,xP(O:~) ,XP~ (0:6) , $ xP2(0:6) ,xP3(0:6> ,~~4(0:6) ,YP(O:6) ,YP~ (0~6) , $ yP2(0:6) ,~~3(0:6) ,~~4(0:6) ,xPP(O:~) ,XPP~ (0~6) , $ xpP2(0:6) ,XPP~ (0:6) ,xpp4(0:6) ,YPP (0:6) ,YPP~ (0:6) , $ YPP2 (0:6) , YPP3 (0: 6) , YPP4 (0: 6) , CD (5) , CL (5) 

COMMON/BLOK~/ TX,TXP,TY,TYP,K,Q,H,L,J,I 

REAL K 
DIMENSION TX (3) , TXP (3) , TY (3) , TYP (3) , K (4) , Q (4) 

C 

DIFF=DELO+TX (2) -TX (3) 

XMAS= PIDEN* P PI/^.) * (DIFF**~. ) +DIFF* (ALo+O. S*DELO+TE 

$ +O.~*(~.-PI)*(TX(~)-TX(~)))+~.~*(~.-PI)*((~.~*DIFF 

$ +TE) **2. ) ) +EM 

YMAS=XMAS 

C 

CD(l)=-0.5 

CL(1)-0.2 

DO 10 KOUNT=2,N 

CD (KOUNT) ~0.085 

CL (KOUNT) =O . 08 

10 CONTINUE 

C 

TERM1-UCT+TYP (2) -TYP (3) 

TERM2=UC+TXP (2) -TXP (3) 

C 

XPP1 (J) =FCNl (XMAS, TERM1, TERM2 , TX (1) , TX (2) , TX (3) , TXP (1) , $ TXP (2) , TXP (3) , TYP (2) , TYP (3) , CD (J) , CL ( J) , FDX (J) , FDY (J) ) 

K (L) =H*XPP 1 ( J) 

C 

YPP1 (J) eFCN2 (YMAS y TERM1 y TERM2 TY (1) STY (2) TY (3) TYP (1) 

TYP (2) , TYP (3) , TXP (2) , TXP (3) , CD (J) JL ( J) , ( J) ) 

RETURN 

END 

FUNCTION FCN~(XMAS,TEM~,TERM~,X~,X~,X~,VX~,VX~,VX~,VY~, 

VY~,CD,CL,DX,DY) 

COMMON/BLOKZ/ XK,YK,FDX,FDY,DENSA 

DIMENS ION FDX (9) , FDY (9) 

RETURN 

END 

FUNCTION FCN~(YMAS,TERM~,TEM~,Y~,Y~,Y~,VY~,VY~,VY~,VX~, 

vX~,CD,CL,DX,DY) 

COMMON/BLOK~/ XK,YK,FDX,FDY,DENSA 

DIMENS ION FDX (9) , FDY (9) 

FCN2= (DENSA* (SQRT ( TERM^ **2. ) + TERM^**^. ) ) ) * (cL*TERM~ 

+CD*TERM~) +DY*SIGN (I., VY~) +YK* (Y 1-2. *~2+~3)) /YMAS 

RETURN 

END 

C 

SUBROUTINE HAMM 

C 

C IMPLEMENTING HAMMINGS METHOD FOR SOLVING GOVERNING EQUATIONS. 

- 71 - 
COMMON/BLOK~/ PI,DELO,PIDEN,ALO,TE,EM,UC,UCT,N 

COMMON/BLOK~/ XK,YK,FDX,FDY,DENSA 

DIMENSION FDX (9) , FDY (9) 

COMMON/BLOK~/ X,Y 

COMMON/BLOK~~/ XP,XP~,XP~,XP~,XP~,YP,YP~,YP~,YP~,YP~, 

~PP,XPP~,XPP~,XPP~,XPP~,YPP,YPP~,YPP~,YPP~,YPP~, 

CD, CL 

DIMENSION ~(0:6,0:2048) ,~(0:6,0:2048) ,xP(O:~) ,XP1(0:6), 

xp2(0:6> ,Xp3(0:6) ,XP4(0:6) ,YP(0:6) ,YP1(0:6), 

Yp2(0:6) ,YP~ (0:6) ,Y~4(0:6) ,XPP (04) ,XPP1(0:6), 

xpp2(0:6) ,xPP3(0:6> ,xPP4(0:6) ,YPP(O:6) ,YPPl(O:6), 

YPP~ (0:6) , YPP3 (0:6> , YPP4 (0:6) , CD(5) ,CL (5) 

DIMENSION PX(O:~) ,PX~ (0:6) ,PXP(O:~) ,PXP~ (0:6) ,PY (0261, 

PY~ (0:6) ,PYP(O:~) ,PYp1(0:6) ,cX(0:6) ,CX~ (0~6) ,cxP(O:~), 

CXP~ (0:6) ,CY (0:6) ,CY~ (0:6) ,CYP(O:6) ,CYP1(0:6) ,lfX(0:6), 

M~(0:6) ,MY (O:6) ,MYP(~:~) 

REAL MX,MXP,MY,MYP 

TIME STEP, H, AND EQUATION CONSTANTS. 

INITIAL VALUE LOOP. 

W 3, JzO,N+1 

1.53 

PX(J)=X(J,I) 

CX(J)=X(J, I) 

PXP (J) =XP 1 (J) 

CXP (J) =XP 1 (J) 

PY (J)=Y (J,I) 

CY (J)=Y (J, I) 

PYP (J) =YP1 (J) 

CYP (J) =YP 1 ( J) 

CONTINUE 

BOUNDARY VALUE LOOP. 

DO 9, J=O,N+l,N+l 

MX(J)=O.O 

MXP(J)=O.O 

MY(J)=O.O 

MYP ( J) PO. 0 

CONTINUE 

*** HAMMINGS METHOD *** 

REESTABLISH INITIAL CONDITIONS. 
DO 10, J=O,N+l 

X(J,I)=X(J,IND) 

Y(J,I)=Y(J,IND) 

10 CONTINUE 

12 CONTINUE 

C 

14 DO 50, 103,2047 

DO 20, J=1,N 

C 

C PREDICTOR 

C 

PXP1 (J) =XP4 (J) +A* (2. *X.?Pl (J? -XPP2 (J) +2. *XPP3 (J) ) 

PX1 (J)=X(J, I-3)+~*(2.*~~1 (J)-XP2 (J)+~.*XP~(J)) 

PYP1 (J)=YP4 (J) +A* (2. *YPPl (J) -YPP2 (J) +2. *YPP3 (J)) 

PY1 (J>=Y (J, 1-31 +A* (2. *Ypl (J) -YP2 (J) +2. *YP~ (J)) 

C 

C MODIFIER 

C 

MXP ( J) =PXP 1 ( J) -B* (PXP (J) -CXP (J) ) 

MX ( J) =PX1 (J) -B* (PX (J) -CX ( J) ) 

MYP (J)=PYPl (J) -B*(PYP (J)-CYP (J)) 

MY (J) =PY 1 ( J) -B* (PY (J) -CY (J) ) 

20 CONTINUE 

C 

C FUNCTION EVALUATIONS 

C 

DO 30, J=l,N 

DIFF=DELO+MX (J) -MX (J-1) 

XMAS-PIDEN* ( (P1/8.) * (DIFF**~. ) +DIFF* (ALo+O. ~*DELo+~. 5* (1. 

$ -PI>*(MX(J)-M(J-1)))+0.5*(4.-PI)*((~.~*DIFF 

$ +TE) **2. ) ) +EM 

YMAS-XMAS 

TERM15UCT+MYP (J) -MYP (J- 1) 

TERM2=UC+MXP ( J) -MXP (J- 1) 

C 

XPP ( J) = (DENSA* (SQRT ( (TERM1 **2. ) + TERM^""^. ) ) ) * (CL (J) 

$ *TERM~+CD(J) TERM^) +mx (J) *SIGN (1. ,MXP (J)) 

$ +xK*(MX(J+~)-~.*MX(J)+MX(J-~)))/ 

$ XMAS 

C 

YPP ( J) = (DENSA* (SQRT ( (TERMlk*2. ) + (TERM2**2. ) ) ) * (CL (J) 

$ *TERM2+CD(J) TERM^) +FDY (J) *SIGN (1. ,MYP (J) ) 

$ +YK* (MY (~+1) MY MY (J)+MY (J-1))) / 

$ YMAS 

30 CONTINUE 

C 

DO 40, J=l,N 

C 

C CORRECTOR 

C 

CXP1 (J)=C*(g.*x~l (J)-XP3(J)+D*(xPP(J)+2.*XPPl (J) 

$ -XPPZ(J))) 

CX1 (J)=C*(g.*X(J,I>-X(J,I-2)+Dk(CXP1 (J)+~.*xP~ (J) 

$ -XP2(J))) 

CYP1 (J)=C*(~.*YP~ (J)-YP~ (J)+D*(YPP (J)+2.*YPPl (J) 

$ -YPPZ(J))) 

CY1 (J)=C*(g.*Y (J, I)-Y (J, I-~)+D*(CYP~ (J)+~.*YP~ (J) 

$ -YP2(J))) 

C 

C FINAL VALUE - 73 - 
C 

XP (J)=CXPl (J)+AB*(PXPl (J)-CXP1 (J)) 

X(J,1+1)=CXl (J)+AB*(PXl (J)-CXl (J)) 

YP(J)=CYPl (J)+AB*(PYP~ (J)-CYP1 (J)) 

Y (J, 1+1)=cyi (J)+AB*(PYI (J)-cyi (J)) 

40 CONTINUE 

C 

DO 45, J=O,N+l 

PX (J) =PXl (J) 

PXP (J) =PXP 1 (J) 

PY (J) =PY 1 (J) 

PYP (J) =PYPl (J) 

CX (J) =CXl (J) 

CXP (J) =CXP 1 (J) 

CY (J)=CYl (J) 

CYP (J) =CYPl (J) 

XP4 ( J) =XP3 (J) 

XP3 (J) =XP2 (J) 

XP2(J)=XP1 (J) 

XP 1 (J) =XP (J) 

YP4 ( J) =YP3 ( J) 

YP3 (J) =YP2 (J) 

YP2(J)=YP1 (J) 

YP 1 ( J) =YP (J) 

XPP4 (J) =XPP3 (J) 

XPP3 (J) =XPP2 (J) 

XPP2 (J) =XPPl (J) 

XPP 1 (J) =XPP (J) 

YPP4 (J) =YPP3 (J) 

YPP3 (J) =YPP2 (J) 

YPP2 (J) =YPPl (J) 

YPP1 (J) -YPP (J) 

45 CONTINUE 

50 CONTINUE 

60 CONTINUE 

C 

RETURN 

END 

C 

C 

C 

SUBROUTINE STRAIN 

C 

C THIS SUBROUTINE COMPUTES THE STRAIN AMPLITUDE FOR 

C EACH OF THE CONVOLUTE TIPS. 

C 

COMMON/BLOK~/ PI,DELO,PIDEN,ALO,TE,F.M,UC,UCT,N 

COMMON/BLOK3/ X,Y 

COMMON/BLOKS/ Al,E,BOTX2,BOTYZ 

COMKON/BLOK6/ SAMP 

DIMENSION x (0:6,0:2048) , Y (0:6,0: 2048) , SAMP (5) , $ EXMAX (5) , EXMIN (5) , EYMAX (5) , EYMIN (5) , $ SMAX (5) , SMIN (5) 

C 

IF (N .GT. 1) GO TO 2 

AMULT=O .1 

GO TO 4 

2 AMULT=O. 13+0.29*FLOAT (N-2) 

4 BOTX~=AMULT*BOTXP 

BOTY~=AMULT*BOTY~ 
TERMX=E*AI*BOTX~ 

TERMY=E*A~*BOTY~ 

C 

DO 10 J=l,N 

EXMAX ( J) =X (J , 0) /TERMX 

EXMIN (J) =EXMAX ( J) 

EYMAX (J) =Y (J, 0) /TERMY 

EYMIN (J) =EYMAX (J) 

10 CONTINUE 

C 

DO 30 1=1,2064 

DO 20 Jsl,N 

EPX=X (J , I) /TERMX 

EPY=Y (J , I) /TERMY 

EXMAX ( J) 5-1 (EXMAX ( J) , EPX) 

EXMIN ( J) PAMINI (EXMIN ( J) , EPX) 

EYMAX ( J) =AMAX1 (EYMAX (J) , EPY) 

EYMIN (J) =AMIN1 (EYMIN (J) , EPY) 

20 CONTINUE 

30 CONTINUE 

C 

DO 35 J=l,N 

SMAX (J) =EXMAX (J) +EYMAX (J) 

SMIN (J) =EXMIN (J) +EYMIN (J) 

35 CONTINUE 

C 

DO 40 J=l,N 

SAMP (J) = (SMAX ( J) -SMIN (J) ) /2. 

40 CONTINUE 

C 

RETURN 

END 

C 

C 

C 

SUBROUTINE FFT 

C 

C THIS SUBROUTINE TRANSFORMS THE SYSTEM RESPONSE 

C PREDICTED BY THE MODEL INTO THE FREQUENCY DOMAIN 

C SO THAT THE DOMINANT FREQUENCIES MAY BE DETERMINED. 

C 

COMMON/BLOKI/ PI,DELO,PIDEN,ALO,TE,EM,UC,UCT,N 

COMKON/BLOK3/ X,Y 

COMMON/BLOK~/ XMAG,YMAG,OMEG 

DIMENSION ~(0:6,0:2048) ,Y (0:6,0:2048) ,~(0:2048), 

$ XMAG (0: 2048) , YMAG (0: 2048) , OMEG (0 : 2048) 

COMPLEX POWER,U,SUMX,SUMY,TERKX,TERMY,XT,YT 

C 

me2048 

H=0.00001 

T=H*NUM 

C 

DO 50 I=O,NUM 

OMEG (I) = (2. *PI*FLOAT (I) ) /T 

50 CONTINUE 

C 

DO 1000 J=l,N 

C 

DO 100 I=0,1 

A-FLOAT (I) 
POWER=CMPLX (0. , (-2. *P I*A) /FLOAT (NUM) ) 

U (I) =CEXP (POWER) 

100 CONTINUE 

C 

DO 200 1=2,NUM 

u (I) =u (1-1) *u (1) 

200 CONTINUE 

C 

DO 400 K=O,NUM-1 

SUMX=CMPLX(O. ,O.) 

SUMY-CMPLX(0. ,O.) 

LO 300 M=O ,NUM-1 

TERMX-X (J , M) *U (MOD (M*K , NUM) ) 

TERMY =Y (J , M) *U (MOD (M*K , NUM) ) 

SUMX=SUMX+TERMX 

SUMY=SUMY+TERMY 

300 CONTINUE 

XT=H*SUMX 

YT=H*SUMY 

XMAG (K) =CABS (XT) 

YMAG (K) =CABS (YT) 

400 CONTINUE 

C 

C 

CALL FINDFRQ (J) 

C 

C 

1000 CONTINUE 

C 

RETURN 

END 

C 

C 

C 

SUBROUTINE FINDFRQ (J) 

C 

C THIS SUBROUTINE SEARCHES THROUGH THE FREQUENCY RESPONSE 

C DATA AND DETERMINES THE THREE MOST DOMINANT FREQUENCIES 

C IN BOTH THE LONGITUDINAL AND TRANSVERSE DIRECTIONS. 

C 

COMMON/BLOKI/ PI,DELO,PIDEN,ALO,TE,EM,UC,UCT,N 

COMMON/BLOK~/ XMAG,YMAG,OMEG 

COMMON/BLOK~/ XFRQ,YFRQ 

DIMENSION XFRQ (3,s) , YFRQ (3,s) ,XMAG (0:2048) , YMAG (0: 2048) , $ OMEG (0 : 2048) , PK (3) 

C 

P-2. *PI 

C 

DO 10 K=1,3 

XFRQ (K, J) 10.0 

YFRQ (K, J) SO. 0 

PK(K)=O. 0 

10 CONTINUE 

C 

DO 30 I=2,60 

IF (XMAG(I) .LT. XMAG(I-1)) GO TO 30 

IF (XMAG(I) .LT. XMAG(I+~)) GO TO 30 

C 

IF (pK(1) . GT. XMAG(I)) GO TO 22 
PK (2) =PK (1) 

~~(11 =XMAG(I) 

XFRQ (3, J) =XFRQ (2, J) 

XFRQ (2, J) =XFRQ (1, J) 

XFRQ (1, J) =omc (I) 

GO TO 30 

IF (PK(2) .GT. XMAG(1)) GO TO 24 

PK (3) =PK (2) 

PK (2) =XMAG (I) 

XFRQ (3, J) =XFRQ (2, J) 

XFRQ (2, J) -0MEG (I) 

GO TO 30 

IF (PK(3) .GT. XMAG(1)) GO TO 30 

PK (3) =XMAG (I) 

XFRQ (3, J) =OMEG (I) 

CONTINUE 

DO 35 K=1,3 

PK(K)=O.O 

CONTINUE 

IF (PK(I) .GT. YMAG(I)) GO TO 36 

PK (3) =PK (2) 

PK (2) =PK (1) 

PK(~) =YMAG(I) 

YFRQ (3, J) =YFRQ (2, J) 

YFRQ(2, J)=YFRQ(l, J) 

YFRQ (1, J) =OMEG (I) 

GO TO 40 

IF (PK(2) .GT. YMAG(1)) GO TO 38 

PK (3) =PK (2) 

- PK (2) =YMAG (I) 

YFRQ (3, J) =Y FRQ (2, J) 

YFRQ (2, J) =OMEG (I) 

GO TO 40 

IF (~~(31 .GT. YMAG(1)) GO TO 40 

PK(~) =YMAG(I) 

YFRQ (3, J) =OMEG (I) 

CONTINUE 

DO 45 K=1,3 

XFRQ (K, J) =XFRQ (K, J) /P 

YFRQ (K, J) =YFRQ (K, J) /P 

CONTINUE 

RETURN 

END 

SUBROUTINE CYCLE 
C THIS SUBROUTINE CALCULATES THE NUMBER OF CYCLES 

C TO FAILURE FOR EACH CONVOLUTION. 

C 

COMMON/BLOK~/ PI,DELO,PIDEN,ALO,TE,EM,UC,UCT,N 

COMMON/BLOKS/ A1, E, BOTX2, BOTY2 

COMMON/BLOK6/ SAMP 

COMMON/BLOK9/ CFAIL 

DIMENSION CFAIL (5) , SAMP (5) 

REAL NF 

C 

B--0.1 

C=-0 .7 

D=(12.*150000.)/E 

DO 30 J-l,N 

NFm10. 

NUM=I 

FNF=D* ( (2. *NF) **B) + ( (2. *NF) **c) -SAMP ( J) /2. 

FNF1-2. *B*D* ( (2. *NF) ** (B-1 .)) +2. *c* ( (2. *NF) ** (C-1.)) 

DELNF=FNF/FNF~ 

IF (ABS (DELNF) .LE. 0.5) GO TO 20 

NF-NF-DELNF 

NUM=NUM+l 

IF (NUM .GT. 50000) GO TO 40 

GO TO 10 

CFAIL-(J) =NF 

CONTINUE 

GO TO 45 

PRINT*, ' CONVERGENCE PROBLEM IN CYCLE.' 

RETURN 

END 

SUBROUTINE DAMAGE 

C 

C THIS SUBROUTINE COMPUTES THE DAMAGE DONE TO 

C EACH CONVOLUTION DURING SAMPLING PERIOD WHICH 

C INDICATES THE CONVOLUTION MOST LIKELY TO 

C FAIL FIRST. 

C 

COMMON/BLOK~/ PI,DELO,PIDEN,ALO,TE,EM~UC,UCT,N 

COMMON/BLOK8/ XFRQ,YFRQ 

COMMON/BLOK~/ CFAIL 

COMMON/BLOK10/ NFAIL,FMAX,FMIN 

DIMENSION XFRQ(3,S) ,YFRQ(3,5) ,CFAIL(5) ,~(5) 

C 

DO 20 J=l,N 

SWO. 0 

DO 10 K=1,3 

TERM-XFRQ (K, J) +Y FRQ (K, J) 

SUM=SUM+TERM 

10 CONTINUE 

~~~=0.02048*SUM 

D (J) =SUM/CFAIL (J) 

20 CONTINUE 

C 

DTEST=O . 0 

DO 30 J=l,N - 78 - 

IF (D(J) .LT. DTEST) GO TO 30 
DTEST-D (J) 

NFAIL= J 

30 CONTINUE 

C 

XFMAX=O . 0 

XFMIN=1000000.0 

YFMAX-0.0 

YFMIN-1000000.0 

W 40 K=1,3 

IF (XFRQ (K, NFAIL) . EQ. 0.0) GO TO 40 

XFMAX--1 (XFRQ (K, NFAIL) , XFMAX) 

XFMIN-AMIN1 (XFRQ (K, NFAIL) , XFMIN) 

40 CONTINUE 

DO 50 K=1,3 

IF (YFRQ(K,NFAIL) .EQ. 0.0) GO TO 50 

Y FMAX-AMAX1 (YFRQ (K , NFAIL) , Y FMAX) 

YFMIN-AMINl (YFRQ (K,NFAIL) , YFMIN) 

50 CONTINUE 

C 

FMAXPAMAXI (XFMAX , Y FMAX) 

FMIN=AMIN~ (XFMIN, YFMIN) 

C 

RETURN 

END 

C 

C 

C 

SUBROUTINE RESULT 

C 

C THIS SUBROUTINE PRODUCES THE DESIRED OUTPUT. 

C 

COMMON/BLOK~/ PI,DELO,PIDEN,ALO,TE,EM,UC,UCT,N 

COMMON/BLOK9/ CFAIL 

COMMON/BLOKlO/ NFAIL,FMAX,FMIN 

DIMENS ION CFAIL (5) 

C 

TMAX-CFAIL (NFAIL) /FMIN 

TMIN-CFAIL (NFAIL) /FMAX 

C 

PRINT*, '** NUMBER OF CYCLES TO FAILURE:' 

PRINT*, ** , CFAIL (NFAIL) 

PRINT*, I**' 

PRINT*, '** MINIMUM TIME TO FAILURE:' 

PRINT*, ' ** ' , TMIN, ' SEC. ' 

PRINT*, I**' 

PRINT*, I** MAXIMUM TIME TO FAILURE:' 

PRINT*, I** ',TMAX,' SEC.' 

C 

RETURN 

END 
APPENDIX I11 

SAMPLE DATA FROM PROGRAM 

The following charts show the cycles to failure, (Nf), taken from several 

test runs of the computer program. In each case the fluid flowing through the 

bellows is taken to be water with a density of 0.036 1bm/in3. Nc is the 

number of convol utes and Uc i s the flow velocity in ft/s. 

Bellows 1 

Convolute pitch = 0.34 in. 

Convolute height = 0.492 in. 

Ply thickness = 0.015 in. 

Number of plys = 1 

Bellows inside diameter = 4.55 in. 

Elastic modulus = 29 x lo6 psi 

Material density = 0.286 1 bm/in3 

Cycles to Failure 
Be1 1 ows 2 

Convolute pitch = 0.3 in. 

Convolute height = 0.395 in. 

Ply thickness = 0.01 in. 

Number of plys = 2 

Be1 lows inside diameter = 4.47 in. 

Elastic modulus = 29 x lo6 psi 

Material density = 0.286 1 bm/in3 

Cycles to Failure 
MSFC TEST BEW 

FROM REFERENCE 2, 

BELLOWS 

NO* 

5002-1 

5002 -2 

5002 -3 

5002-4 

5002-5 

5002-6 

5005-1 

5005-2 

5 006-1 0 

5009-1 

5011-1 

5013-2 

3013-3 

5028-1 

5028-2 

5034-1 

5034-2 

5034-3 

5034-4 

5034-5 

5034-6 

5 034-7 

5034-8 

5034-10 

5034-13 

5034-15 

AG#l 

AG#2 

AG#3 

5035-1 

5035-2 

5035-3 

Vp 

('Ps) 

62 

72 

48 

56 

70 

70 

33.2 

65.2 

77.7 

40.5 

4 9.4 

50 

45 

65.8 

66.6 

82.5 

84 

84,2 

8444 

65,s 

85.5 

85 

85.8 

70.9 

79.5 

69 

57.8 

74 -8 

80.5 

113 

110 

94 

TIME TO 

FAILURE 

(SEC) 

148-158 

525-535 

313-333 

2255-2309 

10460-10585 

730-736 

1841-1863 

75-100 

1500-1530 

2728-2743 

3107-3133 

2007-2027 

749-760 

787 

150 

1672 

1427-1467 

3391-3402 

3683 

353-373 

226-256 

2199 

560-593 

3580-3760 

1375-1495 

4180 

1604-1634 

20099-20124 

11248-11278 

197 6-1996 

437-477 

7916-8036 

CYCLES TO 

FAILUP 

(no 

.17-a18 

.70-a71 

e34m.36 

2.5-2 6 

16.0-16.1 

950.96 

.78--79 

HI1.7 

1.8 

1.9 

1 3 

r47mr48 

NII 

.I6 

2.3 

1*9-2*0 

S#O 

Se3 

.37-*39 

e28-*32 

2, 9 

.74-.78 

3.9-4.1 

1,6-1.8 

4.5 

1.7 

---- 

---- 

4 .4-4.5 

.98-1.1 

16.3-16.6 

MATERIAL 

321 

321 

321 

321 

321 

321 

321 

321 4 

321 

321 

321 

321 

321 

21-6-9 

21-6-9 

21-6-9 

21-6-9 

21-6-9 

21-6-9 

21-6-9 

21-6-9 

21-6-9 

21-6-9 

21-6-9 

21.6-9 

21-6-9 

21-6-9 

21-6-9 

21-6-9 

21-6-9 

21-6-9 

21-6-9 

VARIATION 

ELBOW 

ELBOW 

ELBOW 

ELBW 

ELBOW 

ELBW 

ELBOW 

eLBW 

ELBOW 

X 

LXI 

'A 

'X 

% 

X 

X 

ACTUAL 

STRESS 

(PSI) 

29,800 

28,300 

29,000 

27,800 

26,500 

28,100 

28,200 

933,000 

27,900 

27,900 

27,600 

28,000 

28,800 

e3,OOO 

37,000 

31,200 

31,200 

31,000 

31,000 

34,060 

34,200 

31,000 

21,900 

31,000 

31,200 

31,000 

31,400 

--ow 

---- 

31,000 

31,800 

31,000 

d 

1