Flow Induced Vibration

FLOW INDUCED VIBRATIONS IN PIPES, A FINITE ELEMENT APPROACH IVAN GRANT Bachelor of Science in Mechanical Engineering Nagpur University Nagpur, India June, 2006 submitted in partial ful? llment of requirements for the degree MASTERS OF SCIENCE IN MECHANICAL ENGINEERING at the CLEVELAND STATE UNIVERSITY May, 2010 This thesis has been approved for the department of MECHANICAL ENGINEERING and the College of Graduate Studies by: Thesis Chairperson, Majid Rashidi, Ph. D. Department & Date Asuquo B. Ebiana, Ph. D. Department & Date Rama S. Gorla, Ph. D. Department & Date ACKNOWLEDGMENTS I would like to thank my advisor Dr. Majid Rashidi and Dr.Paul Bellini, who provided essential support and assistance throughout my graduate career, and also for their guidance which immensely contributed towards the completion of this thesis. This thesis would not have been realized without their support. I would also like to thank Dr. Asuquo. B. Ebiana and Dr. Rama. S. Gorla for being in my thesis committe e. Thanks are also due to my parents,my brother and friends who have encouraged, supported and inspired me. FLOW INDUCED VIBRATIONS IN PIPES, A FINITE ELEMENT APPROACH IVAN GRANT ABSTRACT Flow induced vibrations of pipes with internal ? uid ? ow is studied in this work.Finite Element Analysis methodology is used to determine the critical ? uid velocity that induces the threshold of pipe instability. The partial di? erential equation of motion governing the lateral vibrations of the pipe is employed to develop the sti? ness and inertia matrices corresponding to two of the terms of the equations of motion. The Equation of motion further includes a mixed-derivative term that was treated as a source for a dissipative function. The corresponding matrix with this dissipative function was developed and recognized as the potentially destabilizing factor for the lateral vibrations of the ? id carrying pipe. Two types of boundary conditions, namely simply-supported and cantilevered were consi dered for the pipe. The appropriate mass, sti? ness, and dissipative matrices were developed at an elemental level for the ? uid carrying pipe. These matrices were then assembled to form the overall mass, sti? ness, and dissipative matrices of the entire system. Employing the ? nite element model developed in this work two series of parametric studies were conducted. First, a pipe with a constant wall thickness of 1 mm was analyzed. Then, the parametric studies were extended to a pipe with variable wall thickness.In this case, the wall thickness of the pipe was modeled to taper down from 2. 54 mm to 0. 01 mm. This study shows that the critical velocity of a pipe carrying ? uid can be increased by a factor of six as the result of tapering the wall thickness. iv TABLE OF CONTENTS ABSTRACT LIST OF FIGURES LIST OF TABLES I INTRODUCTION 1. 1 1. 2 1. 3 1. 4 II Overview of Internal Flow Induced Vibrations in Pipes . . . . . . Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . Objective . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Composition of Thesis . . . . . . . . . . . . . . . . . . . . . . . iv vii ix 1 1 2 2 3 FLOW INDUCED VIBRATIONS IN PIPES, A FINITE ELEMENT APPROACH 2. 1 Mathematical Modelling . . . . . . . . . . . . . . . . . . . . . . . 2. 1. 1 2. 2 Equations of Motion . . . . . . . . . . . . . . . . . . . 4 4 4 12 12 Finite Element Model . . . . . . . . . . . . . . . . . . . . . . . . 2. 2. 1 2. 2. 2 2. 2. 3 Shape Functions . . . . . . . . . . . . . . . . . . . . . Formulating the Sti? ness Matrix for a Pipe Carrying Fluid 14 Forming the Matrix for the Force that conforms the Fluid to the Pipe . . . . . . . . . . . . . . . . . . . . . 21 2. 2. 4 2. 2. 5Dissipation Matrix Formulation for a Pipe carrying Fluid 26 Inertia Matrix Formulation for a Pipe carrying Fluid . 28 III FLOW INDUCED VIBRATIONS IN PIPES, A FINITE ELEMENT APPROACH 31 v 3. 1 Forming Global Sti? ness Matrix from Elemental Sti? ness Matrices . . . . . . . . . . . . . . . . . . . . 31 3. 2 Applying Boundary Conditions to Global Sti? ness Matrix for simply supported pipe with ? uid ? ow . . . . 33 3. 3 Applying Boundary Conditions to Global Sti? ness Matrix for a cantilever pipe with ? uid ? ow . . . . . . . 34 3. 4 MATLAB Programs for Assembling Global Matrices for Simply Supported and Cantilever pipe carrying ? uid . . . . . . . . . . 35 35 36 3. 5 3. 6 MATLAB program for a simply supported pipe carrying ? uid . . MATLAB program for a cantilever pipe carrying ? uid . . . . . . IV FLOW INDUCED VIBRATIONS IN PIPES, A FINITE ELEMENT APPROACH 4. 1 V Parametric Study . . . . . . . . . . . . . . . . . . . . . . . . . . 37 37 FLOW INDUCED VIBRATIONS IN PIPES, A FINITE ELEMENT APPROACH 5. 1 Tapered Pipe Carrying Fluid . . . . . . . . . . . . . . . . . . . . 42 42 47 50 50 51 54 MATLAB program for Simply Supported Pipe Carrying Fluid . . MATLAB Program for Cantilever Pipe Carrying Fluid . . . . . . MATLAB Program for Tapered Pipe Carrying Flu id . . . . . . 54 61 68 VI RESULTS AND DISCUSSIONS 6. 1 6. 2 Contribution of the Thesis . . . . . . . . . . . . . . . . . . . . . Future Scope . . . . . . . . . . . . . . . . . . . . . . . . . . . . . BIBLIOGRAPHY Appendices 0. 1 0. 2 0. 3 vi LIST OF FIGURES 2. 1 2. 2 Pinned-Pinned Pipe Carrying Fluid * . . . . . . . . . . . . . . Pipe Carrying Fluid, Forces and Moments acting on Elements (a) Fluid (b) Pipe ** . . . . . . . . . . . . . . . . . . . . . . . . . 5 5 7 9 10 11 13 14 15 16 17 21 33 34 36 2. 3 2. 4 2. 5 2. 6 2. 7 2. 8 2. 9 Force due to Bending . . . . . . . . . . . . . . . . . . . . . . . . .Force that Conforms Fluid to the Curvature of Pipe . . . . . Coriolis Force . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Inertia Force . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Pipe Carrying Fluid . . . . . . . . . . . . . . . . . . . . . . . . . . Beam Element Model . . . . . . . . . . . . . . . . . . . . . . . . . Relationship between Stress and Stra in, Hooks Law . . . . . . 2. 10 Plain sections remain plane . . . . . . . . . . . . . . . . . . . . . 2. 11 Moment of Inertia for an Element in the Beam . . . . . . . . . 2. 12 Pipe Carrying Fluid Model . . . . . . . . . . . . . . . . . . . . . 3. 1 3. 2 3. 4. 1 Representation of Simply Supported Pipe Carrying Fluid . . Representation of Cantilever Pipe Carrying Fluid . . . . . . . Pinned-Free Pipe Carrying Fluid* . . . . . . . . . . . . . . . . . Reduction of Fundamental Frequency for a Pinned-Pinned Pipe with increasing Flow Velocity . . . . . . . . . . . . . . . . 4. 2 Shape Function Plot for a Cantilever Pipe with increasing Flow Velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4. 3 Reduction of Fundamental Frequency for a Cantilever Pipe with increasing Flow Velocity . . . . . . . . . . . . . . . . . . . . 5. 1 Representation of Tapered Pipe Carrying Fluid . . . . . . . 39 40 41 42 vii 5. 2 6. 1 Introducing a Taper in the Pipe Carrying Fluid . . . . . . . . Representation of Pipe Carrying Fluid and Tapered Pipe Carrying Fluid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 47 viii LIST OF TABLES 4. 1 Reduction of Fundamental Frequency for a Pinned-Pinned Pipe with increasing Flow Velocity . . . . . . . . . . . . . . . . 38 4. 2 Reduction of Fundamental Frequency for a Pinned-Free Pipe with increasing Flow Velocity . . . . . . . . . . . . . . . . . . . . 40 5. 1 Reduction of Fundamental Frequency for a Tapered pipe with increasing Flow Velocity . . . . . . . . . . . . . . . . . . . . . . 46 6. 1 Reduction of Fundamental Frequency for a Tapered Pipe with increasing Flow Velocity . . . . . . . . . . . . . . . . . . . . . . . 48 6. 2 Reduction of Fundamental Frequency for a Pinned-Pinned Pipe with increasing Flow Velocity . . . . . . . . . . . . . . . . 49 ix CHAPTER I INTRODUCTION 1. 1 Overview of Internal Flow Induced Vibrations in Pipes The ? ow of a ? uid through a pipe can impose pressures on the walls of the pipe c ausing it to de? ect under certain ? ow conditions. This de? ection of the pipe may lead to structural instability of the pipe.The fundamental natural frequency of a pipe generally decreases with increasing velocity of ? uid ? ow. There are certain cases where decrease in this natural frequency can be very important, such as very high velocity ? ows through ? exible thin-walled pipes such as those used in feed lines to rocket motors and water turbines. The pipe becomes susceptible to resonance or fatigue failure if its natural frequency falls below certain limits. With large ? uid velocities the pipe may become unstable. The most familiar form of this instability is the whipping of an unrestricted garden hose.The study of dynamic response of a ? uid conveying pipe in conjunction with the transient vibration of ruptured pipes reveals that if a pipe ruptures through its cross section, then a ? exible length of unsupported pipe is left spewing out ? uid and is free to whip about and im pact other structures. In power plant plumbing pipe whip is a possible mode of failure. A 1 2 study of the in? uence of the resulting high velocity ? uid on the static and dynamic characteristics of the pipes is therefore necessary. 1. 2 Literature Review Initial investigations on the bending vibrations of a simply supported pipe containing ? id were carried out by Ashley and Haviland[2]. Subsequently,Housner[3] derived the equations of motion of a ? uid conveying pipe more completely and developed an equation relating the fundamental bending frequency of a simply supported pipe to the velocity of the internal ? ow of the ? uid. He also stated that at certain critical velocity, a statically unstable condition could exist. Long[4] presented an alternate solution to Housner’s[3] equation of motion for the simply supported end conditions and also treated the ? xed-free end conditions. He compared the analysis with experimental results to con? rm the mathematical model.His experi mental results were rather inconclusive since the maximum ? uid velocity available for the test was low and change in bending frequency was very small. Other e? orts to treat this subject were made by Benjamin, Niordson[6] and Ta Li. Other solutions to the equations of motion show that type of instability depends on the end conditions of the pipe carrying ? uid. If the ? ow velocity exceeds the critical velocity pipes supported at both ends bow out and buckle[1]. Straight Cantilever pipes fall into ? ow induced vibrations and vibrate at a large amplitude when ? ow velocity exceeds critical velocity[8-11]. . 3 Objective The objective of this thesis is to implement numerical solutions method, more specifically the Finite Element Analysis (FEA) to obtain solutions for di? erent pipe con? gurations and ? uid ? ow characteristics. The governing dynamic equation describing the induced structural vibrations due to internal ? uid ? ow has been formed and dis- 3 cussed. The governing equatio n of motion is a partial di? erential equation that is fourth order in spatial variable and second order in time. Parametric studies have been performed to examine the in? uence of mass distribution along the length of the pipe carrying ? id. 1. 4 Composition of Thesis This thesis is organized according to the following sequences. The equations of motions are derived in chapter(II)for pinned-pinned and ? xed-pinned pipe carrying ? uid. A ? nite element model is created to solve the equation of motion. Elemental matrices are formed for pinned-pinned and ? xed-pinned pipe carrying ? uid. Chapter(III)consists of MATLAB programs that are used to assemble global matrices for the above cases. Boundary conditions are applied and based on the user de? ned parameters fundamental natural frequency for free vibration is calculated for various pipe con? urations. Parametric studies are carried out in the following chapter and results are obtained and discussed. CHAPTER II FLOW INDUCED VIBRATION S IN PIPES, A FINITE ELEMENT APPROACH In this chapter,a mathematical model is formed by developing equations of a straight ? uid conveying pipe and these equations are later solved for the natural frequency and onset of instability of a cantilever and pinned-pinned pipe. 2. 1 2. 1. 1 Mathematical Modelling Equations of Motion Consider a pipe of length L, modulus of elasticity E, and its transverse area moment I. A ? uid ? ows through the pipe at pressure p and density ? t a constant velocity v through the internal pipe cross-section of area A. As the ? uid ? ows through the de? ecting pipe it is accelerated, because of the changing curvature of the pipe and the lateral vibration of the pipeline. The vertical component of ? uid pressure applied to the ? uid element and the pressure force F per unit length applied on the ? uid element by the tube walls oppose these accelerations. Referring to ? gures (2. 1) and 4 5 Figure 2. 1: Pinned-Pinned Pipe Carrying Fluid * (2. 2),balancing the forces in the Y direction on the ? uid element for small deformations, gives F ? A ? ? ? 2Y = ? A( + v )2 Y ? x2 ? t ? x (2. 1) The pressure gradient in the ? uid along the length of the pipe is opposed by the shear stress of the ? uid friction against the tube walls. The sum of the forces parallel Figure 2. 2: Pipe Carrying Fluid, Forces and Moments acting on Elements (a) Fluid (b) Pipe ** to the pipe axis for a constant ? ow velocity gives 0 0 * Flow Induced Vibrations,Robert D. Blevins,Krieger. 1977,P 289 ** Flow Induced Vibrations,Robert D. Blevins,Krieger. 1977,P 289 6 A ?p + ? S = 0 ? x (2. 2) Where S is the inner perimeter of the pipe, and ? s the shear stress on the internal surface of the pipe. The equations of motions of the pipe element are derived as follows. ?T ? 2Y + ? S ? Q 2 = 0 ? x ? x (2. 3) Where Q is the transverse shear force in the pipe and T is the longitudinal tension in the pipe. The forces on the element of the pipe normal to the pipe axis accelerate the pi pe element in the Y direction. For small deformations, ? 2Y ? 2Y ? Q +T 2 ? F =m 2 ? x ? x ? t (2. 4) Where m is the mass per unit length of the empty pipe. The bending moment M in the pipe, the transverse shear force Q and the pipe deformation are related by ? 3Y ?M = EI 3 ? x ? x Q=? (2. 5) Combining all the above equations and eliminating Q and F yields: EI ? 4Y ? 2Y ? ? ? Y + (? A ? T ) 2 + ? A( + v )2 Y + m 2 = 0 4 ? x ? x ? t ? x ? t (2. 6) The shear stress may be eliminated from equation 2. 2 and 2. 3 to give ? (? A ? T ) =0 ? x (2. 7) At the pipe end where x=L, the tension in the pipe is zero and the ? uid pressure is equal to ambient pressure. Thus p=T=0 at x=L, ? A ? T = 0 (2. 8) 7 The equation of motion for a free vibration of a ? uid conveying pipe is found out by substituting ? A ? T = 0 from equation 2. 8 in equation 2. 6 and is given by the equation 2. EI ? 2Y ? 2Y ? 4Y ? 2Y +M 2 =0 + ? Av 2 2 + 2? Av ? x4 ? x ? x? t ? t (2. 9) where the mass per unit length of the pi pe and the ? uid in the pipe is given by M = m + ? A. The next section describes the forces acting on the pipe carrying ? uid for each of the components of eq(2. 9) Y F1 X Z EI ? 4Y ? x4 Figure 2. 3: Force due to Bending Representation of the First Term in the Equation of Motion for a Pipe Carrying Fluid 8 The term EI ? Y is a force component acting on the pipe as a result of bending of ? x4 the pipe. Fig(2. 3) shows a schematic view of this force F1. 4 9 Y F2 X Z ?Av 2 ? 2Y ? x2 Figure 2. : Force that Conforms Fluid to the Curvature of Pipe Representation of the Second Term in the Equation of Motion for a Pipe Carrying Fluid The term ? Av 2 ? Y is a force component acting on the pipe as a result of ? ow ? x2 around a curved pipe. In other words the momentum of the ? uid is changed leading to a force component F2 shown schematically in Fig(2. 4) as a result of the curvature in the pipe. 2 10 Y F3 X Z 2? Av ? 2Y ? x? t Figure 2. 5: Coriolis Force Representation of the Third Term in t he Equation of Motion for a Pipe Carrying Fluid ? Y The term 2? Av ? x? t is the force required to rotate the ? id element as each point 2 in the span rotates with angular velocity. This force is a result of Coriolis E? ect. Fig(2. 5) shows a schematic view of this force F3. 11 Y F4 X Z M ? 2Y ? t2 Figure 2. 6: Inertia Force Representation of the Fourth Term in the Equation of Motion for a Pipe Carrying Fluid The term M ? Y is a force component acting on the pipe as a result of Inertia ? t2 of the pipe and the ? uid ? owing through it. Fig(2. 6) shows a schematic view of this force F4. 2 12 2. 2 Finite Element Model Consider a pipeline span that has a transverse de? ection Y(x,t) from its equillibrium position.The length of the pipe is L,modulus of elasticity of the pipe is E,and the area moment of inertia is I. The density of the ? uid ? owing through the pipe is ? at pressure p and constant velocity v,through the internal pipe cross section having area A. Flow of the ? uid through the de? ecting pipe is accelerated due to the changing curvature of the pipe and the lateral vibration of the pipeline. From the previous section we have the equation of motion for free vibration of a ? uid convering pipe: EI ? 2Y ? 2Y ? 2Y ? 4Y + ? Av 2 2 + 2? Av +M 2 =0 ? x4 ? x ? x? t ? t (2. 10) 2. 2. 1 Shape Functions The essence of the ? ite element method,is to approximate the unknown by an expression given as n w= i=1 Ni ai where Ni are the interpolating shape functions prescribed in terms of linear independent functions and ai are a set of unknown parameters. We shall now derive the shape functions for a pipe element. 13 Y R R x L2 L L1 X Figure 2. 7: Pipe Carrying Fluid Consider an pipe of length L and let at point R be at distance x from the left end. L2=x/L and L1=1-x/L. Forming Shape Functions N 1 = L12 (3 ? 2L1) N 2 = L12 L2L N 3 = L22 (3 ? 2L2) N 4 = ? L1L22 L Substituting the values of L1 and L2 we get (2. 11) (2. 12) (2. 13) (2. 14) N 1 = (1 ? /l)2 (1 + 2x/l) N 2 = (1 ? x/l)2 x/l N 3 = (x/l)2 (3 ? 2x/l) N 4 = ? (1 ? x/l)(x/l)2 (2. 15) (2. 16) (2. 17) (2. 18) 14 2. 2. 2 Formulating the Sti? ness Matrix for a Pipe Carrying Fluid ?1 ?2 W1 W2 Figure 2. 8: Beam Element Model For a two dimensional beam element, the displacement matrix in terms of shape functions can be expressed as ? ? w1 ? ? ? ? ? ?1 ? ? ? [W (x)] = N 1 N 2 N 3 N 4 ? ? ? ? ? w2? ? ? ?2 (2. 19) where N1, N2, N3 and N4 are the displacement shape functions for the two dimensional beam element as stated in equations (2. 15) to (2. 18). The displacements and rotations at end 1 is given by w1, ? and at end 2 is given by w2 , ? 2. Consider the point R inside the beam element of length L as shown in ? gure(2. 7) Let the internal strain energy at point R is given by UR . The internal strain energy at point R can be expressed as: 1 UR = ? 2 where ? is the stress and is the strain at the point R. (2. 20) 15 ? E 1 ? Figure 2. 9: Relationship between Stress and Strain, Hooks Law Also; ? =E Rel ation between stress and strain for elastic material, Hooks Law Substituting the value of ? from equation(2. 21) into equation(2. 20) yields 1 UR = E 2 (2. 21) 2 (2. 22) 16 A1 z B1 w A z B u x Figure 2. 0: Plain sections remain plane Assuming plane sections remain same, = du dx (2. 23) (2. 24) (2. 25) u=z dw dx d2 w =z 2 dx To obtain the internal energy for the whole beam we integrate the internal strain energy at point R over the volume. The internal strain energy for the entire beam is given as: UR dv = U vol (2. 26) Substituting the value of from equation(2. 25) into (2. 26) yields U= vol 1 2 E dv 2 (2. 27) Volume can be expressed as a product of area and length. dv = dA. dx (2. 28) 17 based on the above equation we now integrate equation (2. 27) over the area and over the length. L U= 0 A 1 2 E dAdx 2 (2. 29) Substituting the value of rom equation(2. 25) into equation (2. 28) yields L U= 0 A 1 d2 w E(z 2 )2 dAdx 2 dx (2. 30) Moment of Inertia I for the beam element is given as = dA z Figure 2. 11: Moment of Inertia for an Element in the Beam I= z 2 dA (2. 31) Substituting the value of I from equation(2. 31) into equation(2. 30) yields L U = EI 0 1 d2 w 2 ( ) dx 2 dx2 (2. 32) The above equation for total internal strain energy can be rewritten as L U = EI 0 1 d2 w d2 w ( )( )dx 2 dx2 dx2 (2. 33) 18 The potential energy of the beam is nothing but the total internal strain energy. Therefore, L ? = EI 0 1 d2 w d2 w ( )( )dx 2 dx2 dx2 (2. 34)If A and B are two matrices then applying matrix property of the transpose, yields (AB)T = B T AT (2. 35) We can express the Potential Energy expressed in equation(2. 34) in terms of displacement matrix W(x)equation(2. 19) as, 1 ? = EI 2 From equation (2. 19) we have ? ? w1 ? ? ? ? ? ?1 ? ? ? [W ] = N 1 N 2 N 3 N 4 ? ? ? ? ? w2? ? ? ?2 ? ? N1 ? ? ? ? ? N 2? ? ? [W ]T = ? ? w1 ? 1 w2 ? 2 ? ? ? N 3? ? ? N4 L (W )T (W )dx 0 (2. 36) (2. 37) (2. 38) Substituting the values of W and W T from equation(2. 37) and equation(2. 3 8) in equation(2. 36) yields ? N1 ? ? ? N 2 ? w1 ? 1 w2 ? 2 ? ? ? N 3 ? N4 ? ? ? ? ? ? N1 ? ? ? ? ? w1 ? ? ? ? ?1 ? ? ? ? ? dx (2. 39) ? ? ? w2? ? ? ?2 1 ? = EI 2 L 0 N2 N3 N4 19 where N1, N2, N3 and N4 are the displacement shape functions for the two dimensional beam element as stated in equations (2. 15) to (2. 18). The displacements and rotations at end 1 is given by w1, ? 1 and at end 2 is given by w2 , ? 2. 1 ? = EI 2 L 0 (N 1 ) ? ? ? N 2 N 1 ? w1 ? 1 w2 ? 2 ? ? ? N 3 N 1 ? N4 N1 ? 2 N1 N2 (N 2 )2 N3 N2 N4 N2 N1 N3 N2 N3 (N 3 )2 N4 N3 N1 N4 N2 N4 N3 N4 (N 4 )2 ? w1 ? ? ? ? ? 1 ? ? ? ? ? dx ? ? ?w2? ? ? 2 (2. 40) where ? 2 (N 1 ) ? ? L ? N 2 N 1 ? [K] = ? 0 ? N 3 N 1 ? ? N4 N1 N1 N2 (N 2 )2 N3 N2 N4 N2N1 N3 N2 N3 (N 3 ) 2 N1 N4 ? N4 N3 ? ? N2 N4 ? ? ? dx ? N3 N4 ? ? 2 (N 4 ) (2. 41) N 1 = (1 ? x/l)2 (1 + 2x/l) N 2 = (1 ? x/l)2 x/l N 3 = (x/l)2 (3 ? 2x/l) N 4 = ? (1 ? x/l)(x/l)2 (2. 42) (2. 43) (2. 44) (2. 45) The element sti? ness matrix for the beam is obtained by substit uting the values of shape functions from equations (2. 42) to (2. 45) into equation(2. 41) and integrating every element in the matrix in equation(2. 40) over the length L. 20 The Element sti? ness matrix for a beam element; ? ? 12 6l ? 12 6l ? ? ? ? 2 2? 4l ? 6l 2l ? EI ? 6l ? [K e ] = 3 ? ? l 12 ? 6l 12 ? 6l? ? ? ? ? 2 2 6l 2l ? 6l 4l (2. 46) 1 2. 2. 3 Forming the Matrix for the Force that conforms the Fluid to the Pipe A X ? r ? _______________________ x R Y Figure 2. 12: Pipe Carrying Fluid Model B Consider a pipe carrying ? uid and let R be a point at a distance x from a reference plane AB as shown in ? gure(2. 12). Due to the ? ow of the ? uid through the pipe a force is introduced into the pipe causing the pipe to curve. This force conforms the ? uid to the pipe at all times. Let W be the transverse de? ection of the pipe and ? be angle made by the pipe due to the ? uid ? ow with the neutral axis. ? and ? represent the unit vectors along the X i j ? nd Y axis and r and ? rep resent the two unit vectors at point R along the r and ? ? ? axis. At point R,the vectors r and ? can be expressed as ? r = cos + sin ? i j (2. 47) ? ? = ? sin + cos i j Expression for slope at point R is given by; tan? = dW dx (2. 48) (2. 49) 22 Since the pipe undergoes a small de? ection, hence ? is very small. Therefore; tan? = ? ie ? = dW dx (2. 51) (2. 50) The displacement of a point R at a distance x from the reference plane can be expressed as; ? R = W ? + r? j r We di? erentiate the above equation to get velocity of the ? uid at point R ? ? ? j ? r ? R = W ? + r? + rr ? r = vf ? here vf is the velocity of the ? uid ? ow. Also at time t; r ? d? r= ? dt ie r ? d? d? = r= ? d? dt ? Substituting the value of r in equation(2. 53) yields ? ? ? ? j ? r R = W ? + r? + r (2. 57) (2. 56) (2. 55) (2. 53) (2. 54) (2. 52) ? Substituting the value of r and ? from equations(2. 47) and (2. 48) into equation(2. 56) ? yields; ? ? ? ?j ? R = W ? + r[cos + sin + r? [? sin + cos i j] i j] Sin ce ? is small The velocity at point R is expressed as; ? ? ? i ? j R = Rx? + Ry ? (2. 59) (2. 58) 23 ? ? i ? j ? ? R = (r ? r )? + (W + r? + r? )? ? ? The Y component of velocity R cause the pipe carrying ? id to curve. Therefore, (2. 60) 1 ? ? ? ? T = ? f ARy Ry (2. 61) 2 ? ? where T is the kinetic energy at the point R and Ry is the Y component of velocity,? f is the density of the ? uid,A is the area of cross-section of the pipe. ? ? Substituting the value of Ry from equation(2. 60) yields; 1 ? ? ? ? ? ? ? ? ? T = ? f A[W 2 + r2 ? 2 + r2 ? 2 + 2W r? + 2W ? r + 2rr ] 2 (2. 62) Substituting the value of r from equation(2. 54) and selecting the ? rst,second and the ? fourth terms yields; 1 2 ? ? T = ? f A[W 2 + vf ? 2 + 2W vf ? ] 2 (2. 63) Now substituting the value of ? from equation(2. 51) into equation(2. 3) yields; dW 2 dW dW 1 2 dW 2 ) + vf ( ) + 2vf ( )( )] T = ? f A[( 2 dt dx dt dx From the above equation we have these two terms; 1 2 dW 2 ? f Avf ( ) 2 dx 2? f Avf ( dW dW )( ) dt dx (2. 65) (2. 66) (2. 64) The force acting on the pipe due to the ? uid ? ow can be calculated by integrating the expressions in equations (2. 65) and (2. 66) over the length L. 1 2 dW 2 ? f Avf ( ) 2 dx (2. 67) L The expression in equation(2. 67) represents the force that causes the ? uid to conform to the curvature of the pipe. 2? f Avf ( L dW dW )( ) dt dx (2. 68) 24 The expression in equation(2. 68) represents the coriolis force which causes the ? id in the pipe to whip. The equation(2. 67) can be expressed in terms of displacement shape functions derived for the pipe ? =T ? V ? = L 1 2 dW 2 ? f Avf ( ) 2 dx (2. 69) Rearranging the equation; 2 ? = ? f Avf L 1 dW dW ( )( ) 2 dx dx (2. 70) For a pipe element, the displacement matrix in terms of shape functions can be expressed as ? ? w1 ? ? ? ? ? ?1 ? ? ? [W (x)] = N 1 N 2 N 3 N 4 ? ? ? ? ? w2? ? ? ?2 (2. 71) where N1, N2, N3 and N4 are the displacement shape functions pipe element as stated in equations (2. 15) to (2. 18). The displacements and rotations at end 1 is given by w1, ? 1 and at end 2 is given by w2 , ? . Refer to ? gure(2. 8). Substituting the shape functions determined in equations (2. 15) to (2. 18) ? ? N1 ? ? ? ? ? N 2 ? ? ? ? N1 w1 ? 1 w2 ? 2 ? ? ? N3 ? ? ? ? N4 ? ? w1 ? ? ? ? ? ?1 ? ? ? N 4 ? ? dx (2. 72) ? ? ? w2? ? ? ?2 L 2 ? = ? f Avf 0 N2 N3 25 L 2 ? = ? f Avf 0 (N 1 ) ? ? ? N 2 N 1 ? w1 ? 1 w2 ? 2 ? ? ? N 3 N 1 ? N4 N1 ? 2 N1 N2 (N 2 )2 N3 N2 N4 N2 N1 N3 N2 N3 (N 3 )2 N4 N3 N1 N4 N2 N4 N3 N4 (N 4 )2 ? w1 ? ? ? ? ? 1 ? ? ? ? ? dx ? ? ?w2? ? ? 2 (2. 73) where (N 1 ) ? ? L ? N 2 N 1 ? ? 0 ? N 3 N 1 ? ? N4 N1 ? 2 N1 N2 (N 2 )2 N3 N2 N4 N2 N1 N3 N2 N3 (N 3 ) 2 N1 N4 ? 2 [K2 ] = ? f Avf N4 N3 ? N2 N4 ? ? ? dx ? N3 N4 ? ? 2 (N 4 ) (2. 74) The matrix K2 represents the force that conforms the ? uid to the pipe. Substituting the values of shape functions equations(2. 15) to (2. 18) and integrating it over the length gives us the elemental matrix for the ? 36 3 ? 36 ? ? 4 ? 3 ? Av 2 ? 3 ? [K2 ]e = ? 30l 36 ? 3 36 ? ? 3 ? 1 ? 3 above force. ? 3 ? ? ? 1? ? ? ? ? 3? ? 4 (2. 75) 26 2. 2. 4 Dissipation Matrix Formulation for a Pipe carrying Fluid The dissipation matrix represents the force that causes the ? uid in the pipe to whip creating instability in the system. To formulate this matrix we recall equation (2. 4) and (2. 68) The dissipation function is given by; D= L 2? f Avf ( dW dW )( ) dt dx (2. 76) Where L is the length of the pipe element, ? f is the density of the ? uid, A area of cross-section of the pipe, and vf velocity of the ? uid ? ow. Recalling the displacement shape functions mentioned in equations(2. 15) to (2. 18); N 1 = (1 ? x/l)2 (1 + 2x/l) N 2 = (1 ? x/l)2 x/l N 3 = (x/l)2 (3 ? 2x/l) N 4 = ? (1 ? x/l)(x/l)2 (2. 77) (2. 78) (2. 79) (2. 80) The Dissipation Matrix can be expressed in terms of its displacement shape functions as shown in equations(2. 77) to (2. 80). ? ? N1 ? ? ? ? ? N 2 ? L ? ? D = 2? Avf ? N1 N2 N3 N4 w1 ? 1 w2 ? 2 ? ? ? 0 N3 ? ? ? ? N4 (N 1 ) ? ? ? N 2 N 1 ? w1 ? 1 w2 ? 2 ? ? ? N 3 N 1 ? N4 N1 ? 2 ? ? w1 ? ? ? ? ? ?1 ? ? ? ? ? dx ? ? ? w2? ? ? ?2 (2. 81) N1 N2 (N 2 )2 N3 N2 N4 N2 N1 N3 N2 N3 (N 3 )2 N4 N3 N1 N4 N2 N4 N3 N4 (N 4 )2 L 2? f Avf 0 ? w1 ? ? ? ? ? 1 ? ? ? ? ? dx ? ? ?w2? ? ? 2 (2. 82) 27 Substituting the values of shape functions from equations(2. 77) to (2. 80) and integrating over the length L yields; ? ? ? 30 6 30 ? 6 ? ? ? ? 0 6 ? 1? ?Av ? 6 ? ? [D]e = ? ? 30 30 ? 6 30 6 ? ? ? ? ? 6 1 ? 6 0 [D]e represents the elemental dissipation matrix. (2. 83) 28 2. 2. 5Inertia Matrix Formulation for a Pipe carrying Fluid Consider an element in the pipe having an area dA, length x, volume dv and mass dm. The density of the pipe is ? and let W represent the transverse displacement of the pipe. The displacement model for the Assuming the displacement model of the element to be W (x, t) = [N ]we (t) (2. 84) where W is the vector of displacements,[N] is the matrix of shape functions and we is the vecto r of nodal displacements which is assumed to be a function of time. Let the nodal displacement be expressed as; W = weiwt Nodal Velocity can be found by di? erentiating the equation() with time. W = (iw)weiwt (2. 86) (2. 85) Kinetic Energy of a particle can be expressed as a product of mass and the square of velocity 1 T = mv 2 2 (2. 87) Kinetic energy of the element can be found out by integrating equation(2. 87) over the volume. Also,mass can be expressed as the product of density and volume ie dm = ? dv T = v 1 ? 2 ? W dv 2 (2. 88) The volume of the element can be expressed as the product of area and the length. dv = dA. dx (2. 89) Substituting the value of volume dv from equation(2. 89) into equation(2. 88) and integrating over the area and the length yields; T = ? w2 2 ? ?W 2 dA. dx A L (2. 90) 29 ?dA = ?A A (2. 91) Substituting the value of A ?dA in equation(2. 90) yields; Aw2 2 T = ? W 2 dx L (2. 92) Equation(2. 92) can be written as; Aw2 2 T = ? ? W W dx L (2. 93) The Lagr ange equations are given by d dt where L=T ? V (2. 95) ? L ? w ? ? ? L ? w = (0) (2. 94) is called the Lagrangian function, T is the kinetic energy, V is the potential energy, ? W is the nodal displacement and W is the nodal velocity. The kinetic energy of the element †e† can be expressed as Te = Aw2 2 ? ? W T W dx L (2. 96) ? and where ? is the density and W is the vector of velocities of element e. The expression for T using the eq(2. 9)to (2. 21) can be written as; ? ? N1 ? ? ? ? ? N 2? ? ? w1 ? 1 w2 ? 2 ? ? N 1 N 2 N 3 N 4 ? ? ? N 3? ? ? N4 ? ? w1 ? ? ? ? ? ?1 ? ? ? ? ? dx ? ? ? w2? ? ? ?2 Aw2 T = 2 e (2. 97) L 30 Rewriting the above expression we get; ? (N 1)2 ? ? ? N 2N 1 Aw2 ? Te = w1 ? 1 w2 ? 2 ? ? 2 L ? N 3N 1 ? N 4N 1 ? N 1N 2 N 1N 3 N 1N 4 w1 ? ? 2 (N 2) N 2N 3 N 2N 4? ? ? 1 ? ? ? ? ? dx ? N 3N 2 (N 3)2 N 3N 4? ?w2? ? 2 N 4N 2 N 4N 3 (N 4) ? 2 (2. 98) Recalling the shape functions derived in equations(2. 15) to (2. 18) N 1 = (1 ? x/l)2 (1 + 2x/l) N 2 = (1 ? x/l)2 x/l N 3 = (x/l)2 (3 ? 2x/l) N 4 = ? (1 ? x/l)(x/l)2 (2. 9) (2. 100) (2. 101) (2. 102) Substituting the shape functions from eqs(2. 99) to (2. 102) into eqs(2. 98) yields the elemental mass matrix for a pipe. ? ? 156 22l 54 ? 13l ? ? ? ? 2 2? ? 22l 4l 13l ? 3l ? Ml ? [M ]e = ? ? ? 420 ? 54 13l 156 ? 22l? ? ? ? 2 2 ? 13l ? 3l ? 22l 4l (2. 103) CHAPTER III FLOW INDUCED VIBRATIONS IN PIPES, A FINITE ELEMENT APPROACH 3. 1 Forming Global Sti? ness Matrix from Elemental Sti? ness Matrices Inorder to form a Global Matrix,we start with a 6Ãâ€"6 null matrix,with its six degrees of freedom being translation and rotation of each of the nodes. So our Global Sti? ness matrix looks like this: ? 0 ? ?0 ? ? ? ?0 =? ? ? 0 ? ? ? 0 ? ? 0 ? 0? ? 0? ? ? ? 0? ? ? 0? ? ? 0? ? ? 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 KGlobal (3. 1) 31 32 The two 4Ãâ€"4 element sti? ness matrices are: ? ? 12 6l ? 12 6l ? ? ? ? 4l2 ? 6l 2l2 ? EI ? 6l ? ? e [k1 ] = 3 ? ? l 12 ? 6l 12 ? 6l? ? ? ? ? 2 2 6l 2l ? 6l 4l ? 12 6l ? 12 6l ? (3. 2) ? ? ? ? 2 2? 4l ? 6l 2l ? EI ? 6l ? e [k2 ] = 3 ? ? l 12 ? 6l 12 ? 6l? ? ? ? ? 2 2 6l 2l ? 6l 4l (3. 3) We shall now build the global sti? ness matrix by inserting element 1 ? rst into the global sti? ness matrix. 6l ? 12 6l 0 0? ? 12 ? ? ? 6l 4l2 ? 6l 2l2 0 0? ? ? ? ? ? ? 12 ? 6l 12 ? l 0 0? EI ? ? = 3 ? ? l ? 6l 2 2 2l ? 6l 4l 0 0? ? ? ? ? ? 0 0 0 0 0 0? ? ? ? ? 0 0 0 0 0 0 ? ? KGlobal (3. 4) Inserting element 2 into the global sti? ness matrix ? ? 6l ? 12 6l 0 0 ? ? 12 ? ? ? 6l 4l2 ? 6l 2l2 0 0 ? ? ? ? ? ? ? EI 12 ? 6l (12 + 12) (? 6l + 6l) ? 12 6l ? ? KGlobal = 3 ? ? l ? 6l 2 2 2 2? ? 2l (? 6l + 6l) (4l + 4l ) ? 6l 2l ? ? ? ? ? 0 0 ? 12 ? 6l 12 ? 6l? ? ? ? ? 2 2 0 0 6l 2l ? 6l 4l (3. 5) 33 3. 2 Applying Boundary Conditions to Global Sti? ness Matrix for simply supported pipe with ? uid ? ow When the boundary conditions are applied to a simply supported pipe carrying ? uid, the 6Ãâ€"6 Global Sti? ess Matrix formulated in eq(3. 5) is mo di? ed to a 4Ãâ€"4 Global Sti? ness Matrix. It is as follows; Y 1 2 X L Figure 3. 1: Representation of Simply Supported Pipe Carrying Fluid ? ? 4l2 ?6l 2l2 0 KGlobalS ? ? ? ? EI 6l (12 + 12) (? 6l + 6l) 6l ? ? ? = 3 ? ? l ? 2l2 (? 6l + 6l) (4l2 + 4l2 ) 2l2 ? ? ? ? ? 2 2 0 6l 2l 4l (3. 6) Since the pipe is supported at the two ends the pipe does not de? ect causing its two translational degrees of freedom to go to zero. Hence we end up with the Sti? ness Matrix shown in eq(3. 6) 34 3. 3 Applying Boundary Conditions to Global Sti? ness Matrix for a cantilever pipe with ? id ? ow Y E, I 1 2 X L Figure 3. 2: Representation of Cantilever Pipe Carrying Fluid When the boundary conditions are applied to a Cantilever pipe carrying ? uid, the 6Ãâ€"6 Global Sti? ness Matrix formulated in eq(3. 5) is modi? ed to a 4Ãâ€"4 Global Sti? ness Matrix. It is as follows; ? (12 + 12) (? 6l + 6l) ? 12 6l ? KGlobalS ? ? ? ? ?(? 6l + 6l) (4l2 + 4l2 ) ? 6l 2l2 ? EI ? ? = 3 ? ? ? l ? ?12 ? 6l 12 ? 6l? ? ? ? 6l 2l2 ? 6l 4l2 (3. 7) Since the pipe is supported at one end the pipe does not de? ect or rotate at that end causing translational and rotational degrees of freedom at that end to go to zero.Hence we end up with the Sti? ness Matrix shown in eq(3. 8) 35 3. 4 MATLAB Programs for Assembling Global Matrices for Simply Supported and Cantilever pipe carrying ? uid In this section,we implement the method discussed in section(3. 1) to (3. 3) to form global matrices from the developed elemental matrices of a straight ? uid conveying pipe and these assembled matrices are later solved for the natural frequency and onset of instability of a cantlilever and simply supported pipe carrying ? uid utilizing MATLAB Programs. Consider a pipe of length L, modulus of elasticity E has ? uid ? wing with a velocity v through its inner cross-section having an outside diameter od,and thickness t1. The expression for critical velocity and natural frequency of the simply supported pipe carrying ? uid is given by; wn = ((3. 14)2 /L2 ) vc = (3. 14/L) (E ? I/M ) (3. 8) (3. 9) (E ? I/? A) 3. 5 MATLAB program for a simply supported pipe carrying ? uid The number of elements,density,length,modulus of elasticity of the pipe,density and velocity of ? uid ? owing through the pipe and the thickness of the pipe can be de? ned by the user. Refer to Appendix 1 for the complete MATLAB Program. 36 3. 6MATLAB program for a cantilever pipe carrying ? uid Figure 3. 3: Pinned-Free Pipe Carrying Fluid* The number of elements,density,length,modulus of elasticity of the pipe,density and velocity of ? uid ? owing through the pipe and the thickness of the pipe can be de? ned by the user. The expression for critical velocity and natural frequency of the cantilever pipe carrying ? uid is given by; wn = ((1. 875)2 /L2 ) (E ? I/M ) Where, wn = ((an2 )/L2 ) (EI/M )an = 1. 875, 4. 694, 7. 855 vc = (1. 875/L) (E ? I/? A) (3. 11) (3. 10) Refer to Appendix 2 for the complete MATLAB Program. 0 * Flow Induced Vibrat ions,Robert D.Blevins,Krieger. 1977,P 297 CHAPTER IV FLOW INDUCED VIBRATIONS IN PIPES, A FINITE ELEMENT APPROACH 4. 1 Parametric Study Parametric study has been carried out in this chapter. The study is carried out on a single span steel pipe with a 0. 01 m (0. 4 in. ) diameter and a . 0001 m (0. 004 in. ) thick wall. The other parameters are: Density of the pipe ? p (Kg/m3 ) 8000 Density of the ? uid ? f (Kg/m3 ) 1000 Length of the pipe L (m) 2 Number of elements n 10 Modulus Elasticity E (Gpa) 207 of MATLAB program for the simply supported pipe with ? uid ? ow is utilized for these set of parameters with varying ? uid velocity.Results from this study are shown in the form of graphs and tables. The fundamental frequency of vibration and the critical velocity of ? uid for a simply supported pipe 37 38 carrying ? uid are: ? n 21. 8582 rad/sec vc 16. 0553 m/sec Table 4. 1: Reduction of Fundamental Frequency for a Pinned-Pinned Pipe with increasing Flow Velocity Velocity of Fluid(v) Ve locity Ratio(v/vc) 0 2 4 6 8 10 12 14 16. 0553 0 0. 1246 0. 2491 0. 3737 0. 4983 0. 6228 0. 7474 0. 8720 1 Frequency(w) 21. 8806 21. 5619 20. 5830 18. 8644 16. 2206 12. 1602 3. 7349 0. 3935 0 Frequency Ratio(w/wn) 1 0. 9864 0. 9417 0. 8630 0. 7421 0. 5563 0. 709 0. 0180 0 39 Figure 4. 1: Reduction of Fundamental Frequency for a Pinned-Pinned Pipe with increasing Flow Velocity The fundamental frequency of vibration and the critical velocity of ? uid for a Cantilever pipe carrying ? uid are: ? n 7. 7940 rad/sec vc 9. 5872 m/sec 40 Figure 4. 2: Shape Function Plot for a Cantilever Pipe with increasing Flow Velocity Table 4. 2: Reduction of Fundamental Frequency for a Pinned-Free Pipe with increasing Flow Velocity Velocity of Fluid(v) Velocity Ratio(v/vc) 0 2 4 6 8 9 9. 5872 0 0. 2086 0. 4172 0. 6258 0. 8344 0. 9388 1 Frequency(w) 7. 7940 7. 5968 6. 9807 5. 8549 3. 825 1. 9897 0 Frequency Ratio(w/wn) 1 0. 9747 0. 8957 0. 7512 0. 4981 0. 2553 0 41 Figure 4. 3: Reduction of Fundamental Fr equency for a Cantilever Pipe with increasing Flow Velocity CHAPTER V FLOW INDUCED VIBRATIONS IN PIPES, A FINITE ELEMENT APPROACH E, I v L Figure 5. 1: Representation of Tapered Pipe Carrying Fluid 5. 1 Tapered Pipe Carrying Fluid Consider a pipe of length L, modulus of elasticity E. A ? uid ? ows through the pipe at a velocity v and density ? through the internal pipe cross-section. As the ? uid ? ows through the de? ecting pipe it is accelerated, because of the changing curvature 42 43 f the pipe and the lateral vibration of the pipeline. The vertical component of ? uid pressure applied to the ? uid element and the pressure force F per unit length applied on the ? uid element by the tube walls oppose these accelerations. The input parameters are given by the user. Density of the pipe ? p (Kg/m3 ) 8000 Density of the ? uid ? f (Kg/m3 ) 1000 Length of the pipe L (m) 2 Number of elements n 10 Modulus Elasticity E (Gpa) 207 of For these user de? ned values we introduce a taper in the pipe so that the material property and the length of the pipe with the taper or without the taper remain the same.This is done by keeping the inner diameter of the pipe constant and varying the outer diameter. Refer to ? gure (5. 2) The pipe tapers from one end having a thickness x to the other end having a thickness Pipe Carrying Fluid 9. 8mm OD= 10 mm L=2000 mm x mm t =0. 01 mm ID= 9. 8 mm Tapered Pipe Carrying Fluid Figure 5. 2: Introducing a Taper in the Pipe Carrying Fluid of t = 0. 01mm such that the volume of material is equal to the volume of material 44 for a pipe with no taper. The thickness x of the tapered pipe is now calculated: From ? gure(5. 2) we have †¢ Outer Diameter of the pipe with no taper(OD) 10 mm †¢ Inner Diameter of the pipe(ID) 9. mm †¢ Outer Diameter of thick end of the Tapered pipe (OD1 ) †¢ Length of the pipe(L) 2000 mm †¢ Thickness of thin end of the taper(t) 0. 01 mm †¢ Thickness of thick end of the taper x mm Volume of th e pipe without the taper: V1 = Volume of the pipe with the taper: ? ? L ? 2 V2 = [ (OD1 ) + (ID + 2t)2 ] ? [ (ID2 )] 4 4 3 4 (5. 2) ? (OD2 ? ID2 )L 4 (5. 1) Since the volume of material distributed over the length of the two pipes is equal We have, V1 = V2 (5. 3) Substituting the value for V1 and V2 from equations(5. 1) and (5. 2) into equation(5. 3) yields ? ? ? L ? 2 (OD2 ? ID2 )L = [ (OD1 ) + (ID + 2t)2 ] ? (ID2 )] 4 4 4 3 4 The outer diameter for the thick end of the tapered pipe can be expressed as (5. 4) OD1 = ID + 2x (5. 5) 45 Substituting values of outer diameter(OD),inner diameter(ID),length(L) and thickness(t) into equation (5. 6) yields ? 2 ? ? 2000 ? (10 ? 9. 82 )2000 = [ (9. 8 + 2x)2 + (9. 8 + 0. 02)2 ] ? [ (9. 82 )] 4 4 4 3 4 Solving equation (5. 6) yields (5. 6) x = 2. 24mm (5. 7) Substituting the value of thickness x into equation(5. 5) we get the outer diameter OD1 as OD1 = 14. 268mm (5. 8) Thus, the taper in the pipe varies from a outer diameters of 14. 268 mm to 9 . 82 mm. 46The following MATLAB program is utilized to calculate the fundamental natural frequency of vibration for a tapered pipe carrying ? uid. Refer to Appendix 3 for the complete MATLAB program. Results obtained from the program are given in table (5. 1) Table 5. 1: Reduction of Fundamental Frequency for a Tapered pipe with increasing Flow Velocity Velocity of Fluid(v) Velocity Ratio(v/vc) 0 20 40 60 80 100 103. 3487 0 0. 1935 0. 3870 0. 5806 0. 7741 0. 9676 1 Frequency(w) 40. 8228 40. 083 37. 7783 33. 5980 26. 5798 10. 7122 0 Frequency Ratio(w/wn) . 8100 0. 7784 0. 7337 0. 6525 0. 5162 0. 2080 0The fundamental frequency of vibration and the critical velocity of ? uid for a tapered pipe carrying ? uid obtained from the MATLAB program are: ? n 51. 4917 rad/sec vc 103. 3487 m/sec CHAPTER VI RESULTS AND DISCUSSIONS In the present work, we have utilized numerical method techniques to form the basic elemental matrices for the pinned-pinned and pinned-free pipe carrying ? uid. Matlab programs have been developed and utilized to form global matrices from these elemental matrices and fundamental frequency for free vibration has been calculated for various pipe con? gurations and varying ? uid ? ow velocities.Consider a pipe carrying ? uid having the following user de? ned parameters. E, I v L v Figure 6. 1: Representation of Pipe Carrying Fluid and Tapered Pipe Carrying Fluid 47 48 Density of the pipe ? p (Kg/m3 ) 8000 Density of the ? uid ? f (Kg/m3 ) 1000 Length of the pipe L (m) 2 Number of elements n 10 Modulus Elasticity E (Gpa) 207 of Refer to Appendix 1 and Appendix 3 for the complete MATLAB program Parametric study carried out on a pinned-pinned and tapered pipe for the same material of the pipe and subjected to the same conditions reveal that the tapered pipe is more stable than a pinned-pinned pipe.Comparing the following set of tables justi? es the above statement. The fundamental frequency of vibration and the critical velocity of ? uid for a tapered and a pinned-pinned pipe carrying ? uid are: ? nt 51. 4917 rad/sec ? np 21. 8582 rad/sec vct 103. 3487 m/sec vcp 16. 0553 m/sec Table 6. 1: Reduction of Fundamental Frequency for a Tapered Pipe with increasing Flow Velocity Velocity of Fluid(v) Velocity Ratio(v/vc) 0 20 40 60 80 100 103. 3487 0 0. 1935 0. 3870 0. 5806 0. 7741 0. 9676 1 Frequency(w) 40. 8228 40. 083 37. 7783 33. 5980 26. 5798 10. 7122 0 Frequency Ratio(w/wn) 0. 8100 0. 7784 0. 7337 0. 6525 0. 5162 0. 2080 0 9 Table 6. 2: Reduction of Fundamental Frequency for a Pinned-Pinned Pipe with increasing Flow Velocity Velocity of Fluid(v) Velocity Ratio(v/vc) 0 2 4 6 8 10 12 14 16. 0553 0 0. 1246 0. 2491 0. 3737 0. 4983 0. 6228 0. 7474 0. 8720 1 Frequency(w) 21. 8806 21. 5619 20. 5830 18. 8644 16. 2206 12. 1602 3. 7349 0. 3935 0 Frequency Ratio(w/wn) 1 0. 9864 0. 9417 0. 8630 0. 7421 0. 5563 0. 1709 0. 0180 0 The fundamental frequency for vibration and critical velocity for the onset of instability in tapered pipe is approxim ately three times larger than the pinned-pinned pipe,thus making it more stable. 50 6. 1 Contribution of the Thesis Developed Finite Element Model for vibration analysis of a Pipe Carrying Fluid. †¢ Implemented the above developed model to two di? erent pipe con? gurations: Simply Supported and Cantilever Pipe Carrying Fluid. †¢ Developed MATLAB Programs to solve the Finite Element Models. †¢ Determined the e? ect of ? uid velocities and density on the vibrations of a thin walled Simply Supported and Cantilever pipe carrying ? uid. †¢ The critical velocity and natural frequency of vibrations were determined for the above con? gurations. †¢ Study was carried out on a variable wall thickness pipe and the results obtained show that the critical ? id velocity can be increased when the wall thickness is tapered. 6. 2 Future Scope †¢ Turbulence in Two-Phase Fluids In single-phase ? ow,? uctuations are a direct consequence of turbulence developed in ? uid, whe reas the situation is clearly more complex in two-phase ? ow since the ? uctuation of the mixture itself is added to the inherent turbulence of each phase. †¢ Extend the study to a time dependent ? uid velocity ? owing through the pipe. BIBLIOGRAPHY [1] Doods. H. L and H. Runyan †E? ects of High-Velocity Fluid Flow in the Bending Vibrations and Static Divergence of a Simply Supported Pipe†.National Aeronautics and Space Administration Report NASA TN D-2870 June(1965). [2] Ashley,H and G. Haviland †Bending Vibrations of a Pipe Line Containing Flowing Fluid†. J. Appl. Mech. 17,229-232(1950). [3] Housner,G. W †Bending Vibrations of a Pipe Line Containing Flowing Fluid†. J. Appl. Mech. 19,205-208(1952). [4] Long. R. H †Experimental and Theoretical Study of Transverse Vibration of a tube Containing Flowing Fluid†. J. Appl. Mech. 22,65-68(1955). [5] Liu. H. S and C. D. Mote †Dynamic Response of Pipes Transporting Fluids†. J. Eng. for Industry 96,591-596(1974). 6] Niordson,F. I. N †Vibrations of a Cylinderical Tube Containing Flowing Fluid†. Trans. Roy. Inst. Technol. Stockholm 73(1953). [7] Handelman,G. H †A Note on the transverse Vibration of a tube Containing Flowing Fluid†. Quarterly of Applied Mathematics 13,326-329(1955). [8] Nemat-Nassar,S. S. N. Prasad and G. Herrmann †Destabilizing E? ect on VelocityDependent Forces in Nonconservative Systems†. AIAA J. 4,1276-1280(1966). 51 52 [9] Naguleswaran,S and C. J. H. Williams †Lateral Vibrations of a Pipe Conveying a Fluid†. J. Mech. Eng. Sci. 10,228-238(1968). [10] Herrmann. G and R. W.Bungay †On the Stability of Elastic Systems Subjected to Nonconservative Forces†. J. Appl. Mech. 31,435-440(1964). [11] Gregory. R. W and M. P. Paidoussis †Unstable Oscillations of Tubular Cantilevers Conveying Fluid-I Theory†. Proc. Roy. Soc. (London). Ser. A 293,512-527(1966). [12] S. S. Rao †The Finite Element Method in Engineering†. Pergamon Press Inc. 245294(1982). [13] Michael. R. Hatch †Vibration Simulation Using Matlab and Ansys†. Chapman and Hall/CRC 349-361,392(2001). [14] Robert D. Blevins †Flow Induced Vibrations†. Krieger 289,297(1977). Appendices 53 54 0. 1 MATLAB program for Simply Supported Pipe Carrying FluidMATLAB program for Simply Supported Pipe Carrying Fluid. % The f o l l o w i n g MATLAB Program c a l c u l a t e s t h e Fundamental % N a t u r a l f r e q u e n c y o f v i b r a t i o n , f r e q u e n c y r a t i o (w/wn) % and v e l o c i t y r a t i o ( v / vc ) , f o r a % simply supported pipe carrying f l u i d . % I n o r d e r t o perform t h e above t a s k t h e program a s s e m b l e s % E l e m e n t a l S t i f f n e s s , D i s s i p a t i o n , and I n e r t i a m a t r i c e s % t o form G l o b a l M a t r i c e s which are used t o c a l c u l a t e % Fundamental N a t u r a l % Frequency w . lc ; n um elements =input ( ’ Input number o f e l e m e n t s f o r beam : ’ ) ; % num elements = The u s e r e n t e r s t h e number o f e l e m e n t s % i n which t h e p i p e % has t o be d i v i d e d . n=1: num elements +1;% Number o f nodes ( n ) i s e q u a l t o number o f %e l e m e n t s p l u s one n o d e l =1: num elements ; node2 =2: num elements +1; max nodel=max( n o d e l ) ; max node2=max( node2 ) ; max node used=max( [ max nodel max node2 ] ) ; mnu=max node used ; k=zeros (2? mnu ) ;% C r e a t i n g a G l o b a l S t i f f n e s s Matrix o f z e r o s 55 m =zeros (2? nu ) ;% C r e a t i n g G l o b a l Mass Matrix o f z e r o s x=zeros (2? mnu ) ;% C r e a t i n g G l o b a l Matrix o f z e r o s % f o r t h e f o r c e t h a t conforms f l u i d % to the curvature of the % pipe d=zeros (2? mnu ) ;% C r e a t i n g G l o b a l D i s s i p a t i o n Matrix o f z e r o s %( C o r i o l i s Component ) t=num elements ? 2 ; L=2; % T o t a l l e n g t h o f t h e p i p e i n meters l=L/ num elements ; % Length o f an e l e m e n t t1 =. 0001; od = . 0 1 ; i d=od? 2? t 1 % t h i c k n e s s o f t h e p i p e i n meter % outer diameter of the pipe % inner diameter of the pipeI=pi ? ( od? 4? i d ? 4)/64 % moment o f i n e r t i a o f t h e p i p e E=207? 10? 9; roh =8000; rohw =1000; % Modulus o f e l a s t i c i t y o f t h e p i p e % Density of the pipe % d e n s i t y o f water ( FLuid ) M =roh ? pi ? ( od? 2? i d ? 2)/4 + rohw? pi ? . 2 5 ? i d ? 2 ; % mass per u n i t l e n g t h o f % the pipe + f l u i d rohA=rohw? pi ? ( . 2 5 ? i d ? 2 ) ; l=L/ num elements ; v=0 % v e l o c i t y o f t h e f l u i d f l o w i n g t h r o u g h t h e p i p e %v =16. 0553 z=rohA/M i=sqrt ( ? 1); wn= ( ( 3 . 1 4 ) ? 2 /L? 2)? sqrt (E? I /M) % N a t u r a l Frequency vc =(3. 14/L)? sqrt (E?I /rohA ) % C r i t i c a l V e l o c i t y 56 % Assembling G l o b a l S t i f f n e s s , D i s s i p a t i o n and I n e r t i a M a t r i c e s for j =1: nu m elements d o f 1 =2? n o d e l ( j ) ? 1; d o f 2 =2? n o d e l ( j ) ; d o f 3 =2? node2 ( j ) ? 1; d o f 4 =2? node2 ( j ) ; % S t i f f n e s s Matrix Assembly k ( dof1 , d o f 1 )=k ( dof1 , d o f 1 )+ (12? E? I / l ? 3 ) ; k ( dof2 , d o f 1 )=k ( dof2 , d o f 1 )+ (6? E? I / l ? 2 ) ; k ( dof3 , d o f 1 )=k ( dof3 , d o f 1 )+ (? 12? E? I / l ? 3 ) ; k ( dof4 , d o f 1 )=k ( dof4 , d o f 1 )+ (6? E? I / l ? 2 ) ; k ( dof1 , d o f 2 )=k ( dof1 , d o f 2 )+ (6? E?I / l ? 2 ) ; k ( dof2 , d o f 2 )=k ( dof2 , d o f 2 )+ (4? E? I / l ) ; k ( dof3 , d o f 2 )=k ( dof3 , d o f 2 )+ (? 6? E? I / l ? 2 ) ; k ( dof4 , d o f 2 )=k ( dof4 , d o f 2 )+ (2? E? I / l ) ; k ( dof1 , d o f 3 )=k ( dof1 , d o f 3 )+ (? 12? E? I / l ? 3 ) ; k ( dof2 , d o f 3 )=k ( dof2 , d o f 3 )+ (? 6? E? I / l ? 2 ) ; k ( dof3 , d o f 3 )=k ( dof3 , d o f 3 )+ (12? E? I / l ? 3 ) ; k ( dof4 , d o f 3 )=k ( dof4 , d o f 3 )+ (? 6? E? I / l ? 2 ) ; k ( dof1 , d o f 4 )=k ( dof1 , d o f 4 )+ (6? E? I / l ? 2 ) ; k ( dof2 , d o f 4 )=k ( dof2 , d o f 4 )+ (2? E? I / l ) ; k ( dof3 , d o f 4 )=k ( dof3 , d o f 4 )+ (? ? E? I / l ? 2 ) ; k ( dof4 , d o f 4 )=k ( dof4 , d o f 4 )+ (4? E? I / l ) ; % 57 % Matrix a s s e m b l y f o r t h e second term i e % f o r t h e f o r c e t h a t conforms % f l u i d to the curvature of the pipe x ( dof1 , d o f 1 )=x ( dof1 , d o f 1 )+ ( ( 3 6 ? rohA? v ? 2)/30? l ) ; x ( dof2 , d o f 1 )=x ( dof2 , d o f 1 )+ ( ( 3 ? rohA? v ? 2)/30? l ) ; x ( dof3 , d o f 1 )=x ( dof3 , d o f 1 )+ (( ? 36? rohA? v ? 2)/30? l ) ; x ( dof4 , d o f 1 )=x ( dof4 , d o f 1 )+ ( ( 3 ? rohA? v ? 2)/30? l ) ; x ( dof1 , d o f 2 )=x ( dof1 , d o f 2 )+ ( ( 3 ? ohA? v ? 2)/30? l ) ; x ( dof2 , d o f 2 )=x ( dof2 , d o f 2 )+ ( ( 4 ? rohA? v ? 2)/30? l ) ; x ( dof3 , d o f 2 )=x ( dof3 , d o f 2 )+ (( ? 3? rohA? v ? 2)/30? l ) ; x ( dof4 , d o f 2 )=x ( dof4 , d o f 2 )+ (( ? 1? rohA? v ? 2)/30? l ) ; x ( dof1 , d o f 3 )=x ( dof1 , d o f 3 )+ (( ? 36? rohA? v ? 2)/30? l ) ; x ( dof2 , d o f 3 )=x ( dof2 , d o f 3 )+ (( ? 3? rohA? v ? 2)/30? l ) ; x ( dof3 , d o f 3 )=x ( dof3 , d o f 3 )+ ( ( 3 6 ? rohA? v ? 2)/30? l ) ; x ( dof4 , d o f 3 )=x ( dof4 , d o f 3 )+ (( ? 3? rohA? v ? 2)/30? l ) ; x ( dof1 , d o f 4 )=x ( dof1 , d o f 4 )+ ( ( 3 ? rohA? v ? 2)/30? ) ; x ( dof2 , d o f 4 )=x ( dof2 , d o f 4 )+ (( ? 1? rohA? v ? 2)/30? l ) ; x ( dof3 , d o f 4 )=x ( dof3 , d o f 4 )+ (( ? 3? rohA? v ? 2)/30? l ) ; x ( dof4 , d o f 4 )=x ( dof4 , d o f 4 )+ ( ( 4 ? rohA? v ? 2)/30? l ) ; % % D i s s i p a t i o n Matrix Assembly d ( dof1 , d o f 1 )=d ( dof1 , d o f 1 )+ (2? ( ? 30? rohA? v ) / 6 0 ) ; d ( dof2 , d o f 1 )=d ( dof2 , d o f 1 )+ ( 2 ? ( 6 ? rohA? v ) / 6 0 ) ; d ( dof3 , d o f 1 )=d ( dof3 , d o f 1 )+ ( 2 ? ( 3 0 ? rohA? v ) / 6 0 ) ; 58 d ( dof4 , d o f 1 )=d ( dof4 , d o f 1 )+ (2? ( ? 6? rohA? ) / 6 0 ) ; d ( dof1 , d o f 2 )=d ( dof1 , d o f 2 )+ (2? ( ? 6? rohA? v ) / 6 0 ) ; d ( dof2 , d o f 2 )=d ( dof2 , d o f 2 )+ ( 2 ? ( 0 ? ro hA? v ) / 6 0 ) ; d ( dof3 , d o f 2 )=d ( dof3 , d o f 2 )+ ( 2 ? ( 6 ? rohA? v ) / 6 0 ) ; d ( dof4 , d o f 2 )=d ( dof4 , d o f 2 )+ (2? ( ? 1? rohA? v ) / 6 0 ) ; d ( dof1 , d o f 3 )=d ( dof1 , d o f 3 )+ (2? ( ? 30? rohA? v ) / 6 0 ) ; d ( dof2 , d o f 3 )=d ( dof2 , d o f 3 )+ (2? ( ? 6? rohA? v ) / 6 0 ) ; d ( dof3 , d o f 3 )=d ( dof3 , d o f 3 )+ ( 2 ? ( 3 0 ? rohA? v ) / 6 0 ) ; d ( dof4 , d o f 3 )=d ( dof4 , d o f 3 )+ ( 2 ? ( 6 ? rohA? v ) / 6 0 ) ; ( dof1 , d o f 4 )=d ( dof1 , d o f 4 )+ ( 2 ? ( 6 ? rohA? v ) / 6 0 ) ; d ( dof2 , d o f 4 )=d ( dof2 , d o f 4 )+ ( 2 ? ( 1 ? rohA? v ) / 6 0 ) ; d ( dof3 , d o f 4 )=d ( dof3 , d o f 4 )+ (2? ( ? 6? rohA? v ) / 6 0 ) ; d ( dof4 , d o f 4 )=d ( dof4 , d o f 4 )+ ( 2 ? ( 0 ? rohA? v ) / 6 0 ) ; % % I n e r t i a Matrix Assembly m( dof1 , d o f 1 )=m( dof1 , d o f 1 )+ (156? M? l / 4 2 0 ) ; m( dof2 , d o f 1 )=m( dof2 , d o f 1 )+ (22? l ? 2? M/ 4 2 0 ) ; m( dof3 , d o f 1 )=m( dof3 , d o f 1 )+ (54? l ? M/ 4 2 0 ) ; m( d of4 , d o f 1 )=m( dof4 , d o f 1 )+ (? 3? l ? 2? M/ 4 2 0 ) ; m( dof1 , d o f 2 )=m( dof1 , d o f 2 )+ (22? l ? 2? M/ 4 2 0 ) ; m( dof2 , d o f 2 )=m( dof2 , d o f 2 )+ (4? M? l ? 3 / 4 2 0 ) ; m( dof3 , d o f 2 )=m( dof3 , d o f 2 )+ (13? l ? 2? M/ 4 2 0 ) ; m( dof4 , d o f 2 )=m( dof4 , d o f 2 )+ (? 3? M? l ? 3 / 4 2 0 ) ; 59 m( dof1 , d o f 3 )=m( dof1 , d o f 3 )+ (54? M? l / 4 2 0 ) ; m( dof2 , d o f 3 )=m( dof2 , d o f 3 )+ (13? l ? 2? M/ 4 2 0 ) ; m( dof3 , d o f 3 )=m( dof3 , d o f 3 )+ (156? l ? M/ 4 2 0 ) ; m( dof4 , d o f 3 )=m( dof4 , d o f 3 )+ (? 22? l ? 2? M/ 4 2 0 ) ; m( dof1 , d o f 4 )=m( dof1 , d o f 4 )+ (? 13? l ? 2?M/ 4 2 0 ) ; m( dof2 , d o f 4 )=m( dof2 , d o f 4 )+ (? 3? M? l ? 3 / 4 2 0 ) ; m( dof3 , d o f 4 )=m( dof3 , d o f 4 )+ (? 22? l ? 2? M/ 4 2 0 ) ; m( dof4 , d o f 4 )=m( dof4 , d o f 4 )+ (4? M? l ? 3 / 4 2 0 ) ; end k ( 1 : 1 , : ) = [ ] ;% A p p l y i n g Boundary c o n d i t i o n s k(: ,1:1)=[]; k ( ( 2 ? mnu? 2 ) : ( 2 ? mnu? 2 ) , : ) = [ ] ; k ( : , ( 2 ? mnu? 2 ) : ( 2 ? mnu? 2 ) ) = [ ] ; k x(1:1 ,:)=[]; x(: ,1:1)=[]; x ( ( 2 ? mnu? 2 ) : ( 2 ? mnu? 2 ) , : ) = [ ] ; x ( : , ( 2 ? mnu? 2 ) : ( 2 ? mnu? 2 ) ) = [ ] ; x; % G l o b a l Matrix f o r t h e % Force t h a t conforms f l u i d t o p i p e x1=? d(1:1 ,:)=[]; d(: ,1:1)=[]; d ( ( 2 ? mnu? 2 ) : ( 2 ? mnu? 2 ) , : ) = [ ] ; % G l o b a l S t i f f n e s s Matrix 60 d ( : , ( 2 ? mnu? 2 ) : ( 2 ? mnu? 2 ) ) = [ ] ; d d1=(? d ) Kg lobal=k+10? x1 ; m( 1 : 1 , : ) = [ ] ; m( : , 1 : 1 ) = [ ] ; m( ( 2 ? mnu? 2 ) : ( 2 ? mnu? 2 ) , : ) = [ ] ; m( : , ( 2 ? mnu? 2 ) : ( 2 ? mnu? 2 ) ) = [ ] ; m; eye ( t ) ; zeros ( t ) ; H=[? inv (m) ? ( d1 ) ? inv (m)? Kglobal ; eye ( t ) zeros ( t ) ] ; Evalue=eig (H) % E i g e n v a l u e s v r a t i o=v/ vc % V e l o c i t y Ratio % G l o b a l Mass Matrix % G l o b a l D i s s i p a t i o nMatrix i v 2=imag ( Evalue ) ; i v 2 1=min( abs ( i v 2 ) ) ; w1 = ( i v 2 1 ) wn w r a t i o=w1/wn vc % Frequency Ratio % Fundamental N a t u r a l f r e q u e n c y 61 0. 2 MATLAB Program for Cantilever Pipe Carrying Fluid MATLAB Program for Cantilever Pipe Carrying Fluid. % The f o l l o w i n g MATLAB Program c a l c u l a t e s t h e Fundamental % N a t u r a l f r e q u e n c y o f v i b r a t i o n , f r e q u e n c y r a t i o (w/wn) % and v e l o c i t y r a t i o ( v / vc ) , f o r a c a n t i l e v e r p i p e % carrying f l u i d . I n o r d e r t o perform t h e above t a s k t h e program a s s e m b l e s % E l e m e n t a l S t i f f n e s s , D i s s i p a t i o n , and I n e r t i a m a t r i c e s % t o form G l o b a l M a t r i c e s which are used % t o c a l c u l a t e Fundamental N a t u r a l % Frequency w . clc ; num elements =input ( ’ Input number o f e l e m e n t s f o r Pipe : ’ ) ; % num elements = The u s e r e n t e r s t h e number o f e l e m e n t s % i n which t h e p i p e has t o be d i v i d e d . =1: num elements +1;% Number o f nodes ( n ) i s % e q u a l t o num ber o f e l e m e n t s p l u s one n o d e l =1: num elements ; % Parameters used i n t h e l o o p s node2 =2: num elements +1; max nodel=max( n o d e l ) ; max node2=max( node2 ) ; max node used=max( [ max nodel max node2 ] ) ; mnu=max node used ; k=zeros (2? mnu ) ;% C r e a t i n g a G l o b a l S t i f f n e s s Matrix o f z e r o s 62 m =zeros (2? mnu ) ;% C r e a t i n g G l o b a l Mass Matrix o f z e r o s