A Time History Method for Analysing Operational Piping Vibrations

Vibration failure of piping is a serious problem and a matter of concern for safety and reliability of plant operations. Fatigue is the main cause of such failures. Due to the complexity of the phenomenon no closed form design solutions are available. In our study an analytical technique based on the theory of vibrations in the time domain has been presented. Using the inverse theory, the problem has been reduced to a system of Volterra Integral equations to be solved simultaneously at every time step. The solution of the inverse problem may be used in the conventional method to calculate stresses and end reactions which are important from the perspective of engineering design and condition monitoring. The method is robust, simple and can be easily adopted by practicing engineers.


Introduction
Piping witnesses various vibratory loads throughout its life cycle.These vibrations if not controlled will lead to fatigue failures at points of high stress intensity or could even damage the supports.All these could lead to plant outage or even have more severe consequences like fire and loss of human lives [1].Thus it is imperative that the piping system is to be safeguarded against such hazards.
For the standpoint of engineering design adequacy check, dynamic analysis has to be carried out for the piping system for which the forcing function has to be known.This is the conventional method of analysis, which is also termed mathematically as the direct problem [2][3][4].But the major difficulty in the dealing with the vibration problems lies in the estimation of the forcing function.If the exciting forces can be quantified precisely, the system response can be determined with great accuracy by the existing analytical methods.Thus the estimation of the forcing function is essential for carrying out the dynamic analysis and subsequent engineering design check.
Unfortunately this is not readily possible in most cases since the vibrations in an operating pipeline are flow induced.The complexity of flow patterns and the mechanism of force coupling render the determination of the forcing function extremely difficult.In such a scenario data in the form of field vibration measurements in conjunction with some analytical methods can provide a basis for estimating the dynamic force and stress [5][6][7][8][9].The theory of Inverse Problems [2][3][4] invariably forms the basic theoretical framework for such studies.Inverse Theory has found wide applications in the fields of engineering and mathematics.It has become indispensible where the problems are ill-posed in absence of data.In this sense Inverse Theory has got tremendous practical value.
The vibration problems can be studied in both the frequency and the time domains.In the frequency domain the frequency response of the system is studied for the determination of various parameters of interest [7][8][9][10].For a recent work in the frequency domain for a closely related application, one may refer to [8].In the time domain the response time history of different points in the form of observations are taken as the input for the study.
Some good amount of work has been done in the application of inverse problems in the field of hyperbolic partial differential equations or wave equations [5,[11][12][13][14].The determination of point sources from observations has been the main theme of their studies.The present study may be considered as a special case of such applications.However there are some major differences.Concentrated forces in the form of point loads at interior points in the domain have been considered in the previous works, whereas in our case the point sources are the end moments at the boundaries.Also we have considered damping in the system which simulates the real case and thus the treatment is to a great extent different from the previous ones.According to us such study in the time domain has not been done.A few references [7] can be found wherein dynamic stresses have been computed from displacement measurements in the time domain.The method for computing stresses from the displacement measurements has been shown.However the method of force estimation is not provided.Also the topics on existence and uniqueness of the solution which are the key issues for inverse problems have not been addressed.
In our paper we shall present the theory in the time domain and also the numerical scheme for the problem.The numerical algorithm is simple and it can be easily built into any of the common spreadsheet programs with the help of macros.This we believe will find wide application amongst the practicing engineers.

Current Practice-Vibration Screening Criteria
The current practice is the vibration screening criteria method.In this method the vibration response parameters like velocity or displacements are measured in situ and compared against some acceptance criteria.These are in the form of graphs known vibration severity charts [15].For refinery and petrochemical industries, these charts are being extensively used.They are normally found to be conservative.Another widely used criterion is of ASME OM Code [16] a standard followed for nuclear piping.Here the vibration velocity for a piping span between nodes is the criterion.The limiting value of the velocity is determined by the empirical relationship involving coefficients which depend on several parameters like weld arrangements, mass lumping etc.When the peak value of velocity is less that 12.7 mm/sec, it may be assumed that the piping has sufficient dynamic capacity.If the vibration exceeds this level, the guide recommends reviewing the same with more information on the potential reasons of vibrations and improving the vibration levels.
It is seen that all the above methods are conservative and provide a cook book or a go/no-go approach.They only tell us whether the vibrations are within the acceptable levels or not.It is not possible to have a quantitative estimation of the forcing function and the actual stress levels which are essential for a design check.In our work this problem has been studied on the framework of Inverse Theory as mentioned earlier.

Mathematical Background
It has been shown that for a simply supported pipe [vide Figure 1] the response at any location in the span may be determined by the vibration measurements at two distinct points in the span.For the straight span, the excitation source is through moments from the adjoining segments as there are no points of excitation by forces in the span.A distinguishing feature of this method is that no information is required on the natural B.C's.This is remarkable since in the direct formulation, the B.C's govern the solution, whereas in this case they are not playing a role.This is also significant for the fact that practically it is impossible to measure the B.C's.

Notations
In this section we will describe the notations used in the sequel.Please refer to Table 1.

Problem Formulation
The basic configuration is shown in Figure 1 in which the simply supported pipe is excited by moments ( ) ( ) L M t at the ends.The length of the span is L. Considering Bernoulli-Euler formulation and viscous damping, the dynamic equation of motion in the time domain [17,18] is as follows:    Boundary Conditions (B.C's): Equation ( 1) pertains to vibrations without any external loading in the span.It is similar to the free vibration equation.However for our case the excitations are through end moments.This is shown in B.C's (3).In absence of any forcing function in the span the sources of vibrations are through the ends.As mentioned earlier, this is a significant development, since in the earlier studies the point sources of excitation forces have been dealt with.Our study is aimed at the determination of the end moments by observation of the response of some internal points.Then the response at any point in the span can be determined.It also assumed that the system starts from rest (i.e. it has zero initial conditions).
We now express the total displacement function in terms of dynamic and quasi-static components as below.
The quasi-static part can be written in terms of the shape or participation functions ( ) The function Υ is defined as the displacement of the points in the span with a unit positive rotation at end at 0 x = and x L = respectively.Since the system starts from rest we haves ( ) ( ) ( ) ( ) The following can be easily verified.
( ) ( ) ( ) ( ) ( ) ( ) Using ( 7) to (10) we can recast (1) as follows: It is customary to consider damping in terms of dynamic displacements only and hence the last term in ( 13) may be dropped.Equation ( 13) represents a forced vibration problem with a distributed loading for a pipe with clamped ends.Equation ( 12) represents the initial conditions.However this being an inverse problem, the forcing function 0 θ  and L θ  (the rotational accelerations) be- come the unknown quantities which are to be determined.Once they are found out, the problem is transformed into a direct problem and is solvable using commonly used numerical methods.The modal superposition method will be the basis of our study in the sequel.
In line with the modal superposition theory the dynamic displacement may be expressed as the sum of modal components as below: Further we define : ( ) ( ) ( ) ( ) With the above properties we get modal equation from (12) as below: Equation ( 19) is the differential equation for the generalized modal displacement ( ) n q t .This is a second order differential in time variable.Two initial conditions are required for its solution.In our case we have zero displacement and zero velocity at time 0 t = .The solution for ( ) ( ) Here . It is also known as the damped natural frequency.
It is clearly seen that we need to get estimates of the rotational accelerations to obtain ( ) n q t .The modal ac- celeration is obtained differentiating (20) twice and using the below identity: We now define the following terms: ( ) ( ) ( ) , e , , The expression for generalized modal acceleration is The total acceleration which is a sum of dynamic and quasi-static components can be written as Substituting (19) in (20) for ( ) Equation ( 30) is the fundamental equation for our study.It is an integral equation of the second kind [11,19].The right hand side (RHS) quantity represents the acceleration which is the observation.The left hand side (LHS) contains the unknown forcing functions in form of rotational accelerations.Our study will focus on the method of solution for the unknown rotational accelerations.

Solution Method
We will now address the aspects of existence and uniqueness by means of the following propositions.

Proposition 1
For a system as defined by the governing differential Equations ( 13) with B.C's (8), ( 9) and initial conditions (12), the response (i.e.displacement, velocity etc.) at any location x can be obtained from the measurement of acceleration time history at any two interior points.
Proof: We begin with the assumption that ( ) belong to the function space ( ) 0, C T (i.e. the space of continuous functions).
It is shown in Appendix A1 that (29) may be reduced to It is seen that (31) represents a pair of Volterra Integral equations [11,19,20] (one for ( ) and the other for ( ) ).
For the first part of the assertion we need to show that the trivial solution is the only solution.It has been proved in Appendix A2 that both the integral equations have unique fixed points.It is also seen that the trivial solution exists.Hence it is also the unique solution.
As the existence and uniqueness of (31) has been established, we can solve for the individual forcing functions at any time t from the below system: Equations 32 and 33 represent a system of integral equations.The numerical solution method for a single equation may be found in standard texts on numerical methods [20].Thus Equations ( 32) and (33) may be expressed in the matrix form as In an expanded form (34) can be written as The elements of the matrix are as below ( ) ( ) where 1, 2 i = and 0,1 j = and N is the number of modes.The solution t X is obtained from (30) as be- low.
is the RHS vector of the measurements).From the existence of the solution we know that operator ( ) t A is invertible and hence is the unique solution.

Proposition 2
The response obtained at any point is unique and independent of the observation points.This means that if ( ) 1 , v x t is the response calculated on the basis of obser- vations for the set of points ( ) , x x we have 1 2 v v = .Proof: From Proposition 1 we know that the rotational accelerations are determined uniquely.The response is calculated on the basis of the solution of the direct problem (12) which is also unique.The forcing function is the same for all cases.Hence we have 1 2 v v = .

Numerical Method
In this section we shall describe the numerical scheme for the calculation of acceleration forcing function.Let T be the total time interval for our study and T N the number of time steps and N the number of modes.The objective of the scheme is to obtain t X for all the time in- stants t 1 , t 2 ,•••etc.up to T. For convenience will be denoted as i X .The steps are described below: 1) Start with 2) For any r t we define the following quantities ( , The components of the matrix A in (35) is constructed as ( ) The RHS vector U is ( ) (Here i = 1 and 2) 3) Solve for r X from (39).4) Repeat Steps 2 to 3 till T r N =

Determination of Response Variables
We obtain the rotational accelerations as a solution of the inverse problem.Now we can determine the forcing function completely.Thus the problem is transformed into a direct one, which may be solved using existing methods for determining various response quantities like displacement, velocity and stress time histories.For example, Bending Moment, Shear Force and the Bending Stress are calculated as below. , As a measure of structural integrity a mechanical design check against fatigue is required to be carried out using the stress distribution.In the time domain it is customary to apply Rain-Flow Counting Method [21] to determine cumulative usage factor.The value of the usage factor should be less than unity which indicates that the system is safe and no failures from fatigue are expected to occur in its design life.A value of the factor greater than one is an indication of a possibility of failure due to fatigue.However as a crude estimate we may consider maximum zero to peak value of the stress and compare it with the endurance limit.It should be less than the endurance limit to designate a system as safe.
The time histories of velocities and the end reactions can be computed through the direct problem.The end reaction forces should be used for checking mechanical design of the support structure.This will ensure integrity of the pipe supports thereby accounting for an important hazard of a vibrating piping system.
The velocity at a point may be compared against the maximum permissible velocity as per common practices as mentioned earlier.However in view of our detailed analytical method they are not the essential parameters and may be taken as an additional piece of information.

Numerical Simulation and Validation
In order to validate the theory some numerical experiments have been carried out.The problem considered is as follows.
A simply supported pipe is excited through end moments.Two cases have been considered.In case 1 the excitation moment is applied at only one end.In case 2 excitation moments are applied at both ends.For simplicity the harmonic excitation comprising of sine and cosine terms for a few frequencies have been considered for the forcing functions.However any continuous time varying function is permissible.The total time T considered is 200 seconds.The pipe material is steel, size 219 mm outer diameter (O.D), thickness 8.18 mm and the span is 8 m.A fixed damping ratio of 1% has been assumed.Five points numbered 1 to 5 have been defined in the span.Points 1 at x = 0 and 5 at x = L are the boundary points.Points 2, 3 and 4 are interior points at locations 0.25 L, 0.5 L and 0.75 L respectively.These points have been defined for the purpose of specifying the input and output locations.
The direct problem is first solved using the forcing function as the moments using standard software.The dynamic analysis time history module of general purpose Finite Element Analysis (FEA) software has been used for the direct problem.This analysis model will be termed as model D in the sequel.The results of the analysis have been treated as the benchmark.The accelerations from model D have been considered as measurements which are the inputs for our proposed method which is based on Inverse Theory and denoted as model I for reference.Displacements, stresses and end reactions have been considered as the response parameters for comparison with the benchmark.

Results and Discussions
The time step interval has been fixed based on the highest natural frequency.This is done for the purpose of minimizing errors due to integration.For the details on the theory one may refer to standard texts [17,18].Five modes have been considered for the problem.The observation points are 2 and 4 where the acceleration time histories are measured [see Figures 3 and 9].The rotational accelerations are calculated from Equation 37 as per inverse theory.It is seen from Figure 4 that the rotational accelerations are shown at point 1 only.This is due to the fact that in Case 1 the excitations are applied at one end.The other response quantities like end reactions, displacements and stresses are shown in Figures 5- 7.In all cases there is no difference between the results of the two models.In the sequel we shall use the abbreviations TH for time history, ATH for acceleration time history and RTH for rotational time history.
The results for Case 2 are given in Figures 9-13.In this case we have rotational accelerations for both the ends unlike Case 1.Also a very close match between the results of direct and inverse problem is observed similar to Case 1.This is expected since the theoretical solution for the two methods is essentially the same.The difference is basically due to the round off errors.As mentioned earlier, the distinct advantage of the method over the current ones is that quantitative estimate of the stresses and the end reactions are obtained in this method.This is significant from the aspect of condition monitoring and engineering design.The reaction force  estimates will enable us to design the pipe supports, whereas the stresses and displacements will be useful for condition monitoring of the system.

Conclusion
Vibration failure in operational piping is a serious problem and there is a need for a comprehensive study and analysis for its remedial measures.In this sense the proposed study has got a tremendous practical value.A quantitative method with proper mathematical basis has been provided in contrast to the cook book approach.By this method it is possible to quantify stresses, velocities and reaction forces.This gives us a basis for a proper engineering design.The method being simple can be easily adopted by engineers involved in trouble-shooting.Several improvements in the model are in line and planned for future work.These are like inclusion of lumped mass in the span or pipe bends.These will widen the range of application of the method and will be of greater practical use.
Space of continuous functions in [0, T]. u Acceleration measurement time history.0, L K KKernel of the integral equation.TH.Time History.
-function for the th n mode for the clamped pipe.

Figures 2
and show the moment time history for Case 1 and Case 2 respectively.Graph D denotes the input for direct problem model whereas graph I denotes the calculated response for the Inverse Problem.It is seen that the two graphs coincide implying unique correspondence between the Inverse and Direct Problem for our case.
n ωUndamped natural frequency (n th mode).