A Semi-Discretizing Method Based Efficient Model for Fluidelastic Instability Threshold Prediction of Tube Bundles

Fluidelastic instability is destructive in tube bundles subjected to cross flow. Flow channel model proposed by Leaver and Weaver is well used for modeling this problem. However, as the tube motion is supposed to be harmonic, it may not simulate the general dynamic behaviors of tubes. To improve this, a model with arbitrary tube motion is proposed by Hassan and Hayder. While, due to involving in the time delay term, the stability problem cannot be solved by the eigenvalue scheme, and time domain responses of the tube have to be obtained to assess the instability threshold. To overcome this weakness, a new approach based on semi-discretizing method (SDM) is proposed in this study to make the instability threshold be predicted by eigenvalues directly. The motion equation of tube is built with considering the arbitrary tube motion and the time delay between fluid flow and tube vibration. A time delay integral term is derived and the SDM is employed to construct a transfer matrix, which transforms the infinite dimensional eigenvalue problem into a finite one. Hence the stability problem become solvable accordingly. With the proposed method, the instability threshold of a typical square tube array model is predicted, and the influences of system parameters on stability are also discussed. With comparing with prior works, it shows significant efficiency improvement in prediction of the instability threshold of tube bundles.


Introduction
Tube and shell heat exchangers are widely used in conventional and nuclear industries. Tube bundles in heat exchangers are usually exposed to external flows, and the interaction between tube and fluid flow will induce tube vibration. It often leads to damage and even failures of tubes which has become a critical problem of heat exchangers [Shinde (2015)]. The excitation mechanisms of Flow-Induced Vibration (FIV) can be identified as turbulence excitation, vortex shedding, fluidelastic instability and acoustic resonance broadly. As fluidelastic instability may cause severe unacceptable damage and is the most complex and difficult to predict, it has attracted more attentions than other three FIV behaviors. The feature phenomenon of fluidelastic instability is that when the flow velocity exceeds certain threshold, the amplitude of tube vibration rises sharply with a small increment of inlet flow velocity, and it may lead to the tube failure in a short time. Due to the danger of fluidelastic instability of tube bundles, amounts of researches have been carried out to predict and prevent it, and several theoretical models for tube arrays subjected to cross flow were developed. As the fluid forces are difficult to solve theoretically, the early solution was to build a simplified model in terms of extensive experimental data. In Roberts' experiments [Roberts Washington (1966)], transverse fluidelastic instability was observed in the in-flow direction at low flow velocities. Robert made the first attempt at analyzing the phenomenon by the Jet Switch model. Researching on a single row of tube array, Connors [Connors (1970)] proposed quasi-static model. A quantity expression between reduced flow velocity and mass-damping parameters (MDP) is given in the model. Connors did not approach the fluid forces theoretically, but measured them. Takahara (1980, 1981)], Chen [Chen (1983a[Chen ( , 1983b] provided an unsteady model, where the unsteady forces were obtained directly from experiments. Also, there were quasi-steady models developed by Blevins [Blevins (1977)] and Price et al. [Price andPaidoussis (1984, 1986)], where four new constants were to be determined and couldnot account for stiffness-controlled instabilities. A semi-analytical approach to modeling fluidelastic instability was presented by Lever et al. [Lever and Weaver (1982); Weaver (1986a, 1986b)] and Yetisir et al. Weaver (1993a, 1993b)], which was based on damping-controlled mechanism. The model was considered in the case that an elastic tube was confined in a rigid tube bundle. In the original form, the flow channel model of Lever et al. [Lever and Weaver (1982)] idealized the tube as a single degree of freedom system vibrating at its natural frequency, only in the transverse direction, and a time delay was expressed associating with the time taken for the two fluid streams on either side of the tube to readjust to the changing flow-channel configuration as a tube vibrates. Lever et al. Weaver (1986a, 1986b)] modified the flow channel to account for streamwise motion, however, the transverse and streamwise flow motions were analyzed independently of each other. Weaver (1993a, 1993b)] then extended the original model to account for the case of multiple flexible tubes. Since the differential motion equation was linear, the stability problem can be assessed in terms of the roots of the characteristic polynomial. Although Lever and Weaver's model described the system behavior and fluid force theoretically, due to the assumption that tube motion is harmonic and only in the transverse direction, which is not identical to its arbitrary vibration facts, this makes the model not be suitable for nonlinear analysis. Hassan et al. [Hassan and Hayder (2008); Hassan and Weaver (2016); Hassan and Weaver (2017)] modified the semi-analytical model to account for any arbitrary motion, expanding the motion of tube to two dimensions in the plane, and developed a time domain model. Although this model may approach the arbitrary tube free motion with no limitation, it brought in a heavy burden in solving the equations by direct time integration scheme. On account of that only one inlet flow speed case corresponding to a certain MDP can be solved in once numerical computation. Moreover, time domain responses of the tube at different inlet flow speeds have to be obtained for assessing the instability via vibrational amplitudes. It will take a large amount of calculation time and hence limit its application in engineering analysis. To avoid the time domain solution for the Hassan and Hayder's model, the obstacle is that the time delay and unlimited vibration mode leads to an infinite-dimensional eigenvalue problem as the multiplier of characteristic polynomial has no closed form. This paper hence proposes a tube-flow coupling dynamic equation with a constant delay term and a time delay integral term based on Hassan and Hayder's model, which supposes the tube motion is arbitrary. The SDM is then introduced to overcome the challenge via desecrating the time delay integration term. Furthermore, a transform matrix is built to make the stability problem become a finite eigenvalue one. The stability characteristics are presented by the stability map in different MDPs and effects of related parameters on stability are discussed. Finally, the results are compared with the test and theoretical one of prior researches to validate the proposed approach.

Elastic fluid mechanics modeling
In this study, a typical square tube array shown in Fig. 1 is chosen, and the elastic fluid force is approached by reformulated flow channel model proposed by Hassan et al. [Hassan and Hayder (2008)].

Figure 1: Flow channel model
The flow passing tubes is divided into two fluid channels and wake regions as shown in Fig. 1. The flexible tube (tube 5) vibrates in the -plane, while other surrounding tubes are supposed to be rigid and fixed. Since fluid boundaries contact with tubes, it is influenced by motion of tubes as follows. The fluid boundary 1 varies with the motion of the analyzed tube (tube 5). For the fluid boundary 2 , the right end of which varies with the analyzed tube while the other end is fixed. Fluid boundaries 3 , 4 , 5 and 6 are all fixed. The curvilinear coordinate ( ) is measured from the tube center along the curved path of fluid flow and 0 is corresponding to the inlet position of flow channel. For simplifying the modeling, it should be noted that the fluid flow is regarded as incompressible and inviscid, as well as the fluidelastic excitation is assumed to be independent on wake phenomena. The flow is assumed to attach and separate at fixed locations. And these assumptions may not be valid for large amplitude of tube vibration. The conservation of mass equation and momentum equation obtained from the one-dimensional Reynolds equation are in following forms, respectively: where, the fluid flow velocity (s, t), pressure (s, t) and the fluid channel area (s, t) are assumed to vary with position and time respectively, is the one-dimensional curvilinear coordinate along the channel. The fluid boundaries 1 and 2 vary with the vibration of the tube, then the channel area , flow velocity and fluid pressure P fluctuate. The above three parameters can be expressed by the sum formulas of constant components and variable components, which are functions of time and position. The average channel area of the overall fluid channel 0 , the flow velocity 0 and the pressure 0 at the inlet are assumed to be constant [Yetisir and Weaver (1993a)]: Via substituting the Eqs. (3), (4), (5) into (1) and (2), the flow velocity and the pressure distributions along the coordinate can be obtained by solving the equation. By integrating the pressure distribution over the fluid boundary 1 , the force exerting on the analyzed tube (tube 5) can be obtained. Eqs. (4) and (5) are rewritten by utilizing integration by parts as where, ℎ is the fluid resistance coefficient which is related to the arrangement and geometric parameters of tube arrays, supposed not to vary significantly with Reynolds number in the vicinity of the stability threshold for each tube bundle array. As shown in Fig. 2, the pressures in channel 1 and 2 can be calculated from the pressure distributions along the two stream channels, respectively 1 (s, t) and 2 (s, t). The lift and drag force on per unit length along and direction are expressed respectively as and correspond to the positions of flow becoming into attachment and separation, respectively. is the angle between the surface normal and the transverse line at = 0 of the tube. To avoid assessing the stability via time-based solutions, in the following, the fluidelastic force is derived based on Hassan and Hayder's model where SDM can be applied for, and make the instability threshold can be predicted directly in terms of eigenvalues. For small tube responses, that is 0 ≫ ( , ), Eq. (6) becomes And ( , ) is a piecewise function in the curvilinear coordinate, when > , Take the variable substitution ζ = + 2( − ) 0 ⁄ into Eq. (12), where  and (2-14), it is obtained In flow channel model, for any , it gives Combining Eqs. (16) and (8), it yields where 1 , 2 are fluctuating channel area corresponding to flow channel 1 and 2. And 1 , 2 are fluctuating flow velocity corresponding to flow channel 1 and 2. Taking partial derivative with respect to time on both sides of Eq. (17), it is written as follows: For > , it gives The fluid channel area depends on the tube displacement along direction and the time lag. The relationship is given as where, ( ) is the time delay between fluid flow redistribution and the elastic tube motion. So far, there are investigations about time delay in a series of literatures, such as Lever et al. [Lever and Weaver (1982)], Price et al. [Price and Paidoussis (1984)] , Granger [Granger and Paidoussis (1996)], and later researches [Bouzidi, Hassan, Fernandes et al. (2014); El Bouzidi and Hassan (2015); Li and Mureithi (2017)]. In present work, the Lever and Weaver linear model of time delay ( ) is adopted, which is shown as follows, is the relative fluid length coefficient, which related to the phase lag of fluid flow. The phase lag is considered to be closely associated with the fluidelastic instability, which is related to the tube arrangement and pitch ratio. The value of is adopted as 2 in the paper. The artificial time delay function is shown in Fig. 3.
Eqs. (15) and (20) are substituted into Eq. (19) to produce Inserting Eq. (23) into Eq. (8) and integrating it along , the lift and drag force vectors of the tube can be obtained finally. The impacts of downstream wake and phase change due to heat transfer between tubes and fluid flow are not considered, which may lead the model too complicated to solve. The fluid forces are expressed by the state variable of displacement coordinate with the time delay term in this section, then the dynamic equation is established in the following section.

Dynamic equation with time delay of tube bundles
In the stability analysis of tube bundles, it can be considered that only fluidelastic instability force acts [Axisa, Antunes and Villard (1988)]. The Eq. (8) indicates that the fluidelastic force in the direction only depends on the displacement and velocity in the direction, which is considered as an external excitation force. And the stability in the direction of tube bundles only depends on the system structural damping, if the damping is negative, then the system is unstable. Since it is impossible that the system structural damping is negative in practical projects, the tube response along the direction is always stable. The stability problem of tube bundles is then simplified to one dimensional case, i.e., only the stability in direction needs to be solved. As a result, for the stability problem of tube vibration, it is reasonable to assume the tube bundle as an Euler-Bernoulli beam model vibrating in the plane, whose axis is along direction and normal to the plane shown in Fig. 1, and the simplified tube bundle control equation is written as follows: where and are respectively the density and lateral area of the elastic tube, is the structural damping parameter of the system, and are Young's modulus of the beam and the cross-sectional moment of inertia, ( , ) is fluidelastic force in transverse direction. By mode superposition method, extracting preceding N modes, there is ( ) = ∑ ( )

=1
( ). It is demonstrated both by theory and experiments that the critical velocity corresponding to the first order mode is the lowest. Hence intercept the first order mode by Galerkin method and integrate with 1 ( ), then the tube partial differential equation is converted into the ordinary differential equation: 11 + 2 1 1 11 + 1 1 2 1 = ( ) (25) where ( ) is defined as follows, Rearranging Eqs. (23), (25) and (26), the following simplified form is obtained: where = 0 0 ⁄ , ( = 1,2,3,4,5). Making a variable substitution ϑ = 2( − ) 0 ⁄ , we obtain where, Through the above model simplification, the tube dynamic equation is derived, whose right side is composed by a constant delay term and a time delay integral term and the stability problem is transformed to be solvable, where SDM could be applied. The equation is established on the assumption that the tube vibrates in small amplitude in the plane. The new dynamics Eq. (28) takes the first-order mode coordinate as state variable. The details of solving method and process for the proposed model by the SDM are presented in the next session.

Solution of fluidelastic dynamics equation
Autonomous ordinary differential equations (ODEs) have the general form ̇( ) = ( ) ( ) (29) The stability properties are determined by the roots of the characteristic polynomial: if and only if all the characteristic roots have negative real parts, the system is asymptotically stable. Opposite to the characteristic polynomial of autonomous ODEs, for the linear autonomous DDEs in the following form, the characteristic function has infinite number of zeros.
In the DDEs, the sufficient and necessary condition for asymptotic stability is that all the infinite number of characteristic roots have negative real parts. For practical calculations, only approximations can be applied. Stability investigations are often carried out by numerical simulations and the SDM is composed and applied for determining stability criteria for second-order DDEs by Insperger et al. [Insperger and Stépán (2011)]. For the DDE (28), the challenge of fluidelastic instability analysis is that the time delay and unlimited vibration mode lead to an infinite-dimensional eigenvalue problem as the fundamental matrix of characteristic polynomial has no closed form, and it cannot be solved by traditional method. In this section, with the application of SDM, the infinite-dimensional eigenvalue problem is transformed into an approximate finite-dimensional one, and the fluidelastic instability threshold can be predicted by the eigenvalues. SDM has been widely used in solid finite element analysis and computational fluid dynamics, the basic idea of which is to discretize the partial differential equation only along the spatial coordinate without the change of the time coordinate [Insperger and Stépán (2002)]. In this work, the SDM is introduced to construct a transfer matrix, then the stability of the tube-flow coupled control equation with time delay can be determined by its eigenvalues. An overview of stability map will be provided indicating the correlation between the fluidelastic instability threshold and system parameters, which is efficient for predicting of fluidelastic instability. The stability of the dynamic equation is solved by employing SDM as follows. The length of time axis is divided into a series of time intervals, and the period is expressed as = , ∈ . The relationship between and is as Let = [ 1 ̇1] , then Eq. (28) can be rewritten as where, in this work, is an impulse weighting distribution matrix and it is integrable, which satisfies the following equation: At the interval [ , +1 ], (ϑ) = (ϑ). The continuous distribution matrix � (ϑ) can be expressed as a sum of a series of impulse function [Insperger and Stépán (2002)]: Meanwhile, at the interval [ , +1 ], � (ϑ) can also be expressed as a sum of a series of time-dependent distribution [Insperger and Stépán (2002)], With the above approximations, the following approximate form is obtained [Insperger and Stépán (2002)], where − +1 = ( − ( − 1) ), = 0,1, . . . , = 0,1, . . . , Where � i,j = � i,j = 1,2, . . . , = + 1 Assuming is reversible for all , the solution of Eq. (38) is [Insperger and Stépán (2011)], where the constant vector depends on the initial value of the state variable , namely If is irreversible, there is also a corresponding expression, and SDM is still applicable. In our work, it will not be discussed in detail. Eq. (41) is substituted into Eq. (40) to get [Insperger and Stépán (2011)], where the coefficient matrices are Utilizing Eq. (46) we obtain = ⋯ 1 0 0 (48) Then the transfer matrix over the period is given in the form = −1 ⋯ 1 0 (49) By the construction of the transfer matrix , the stability of Eq. (42) then can be assessed by judging whether its eigenvalues are all in the unit circle. If all the eigenvalues are in the unit circle [Lakshmikantham and Trigiante (2002)] , then the solution of the Eq. (32) is stable. Let be the maximum modulus among all eigenvalues for . If > 1, then the system is unstable; if < 1, then the system is asymptotically stable; else when = 1, the system is marginally stable.

Results and discussion
In the view of the damping mechanism for fluidelastic instability of the tube bundle, the system instability occurs when the energy the tube obtained from the flow field is more than that it dissipated. The system damping and flow velocity are important factors affecting system stability. Therefore, the mass-damping parameter (MDP) is introduced as an important factor for fluidelastic instability. The expression of the MDP is where and are respectively the tube mass in per unit length and logarithmic decrement of structural damping, is the tube diameter. The relevant geometrical parameters and material properties in computation are listed in Tab. 1. Fig. 4 shows influences of MDP and reduced flow velocity ( = 0 ⁄ ) on the modulus of the transfer matrix eigenvalues. The stability surface is constructed by computing 50 groups of MDPs respectively in 200 flow velocities and it divides the space into stable and unstable regions. To observe the region where > 1 in ( ) visually, let = 1 as in ( ), then the discrete points corresponding to = 1 compose a curve, namely, the system stability curve. Extracting mass damping parameters and reduced flow velocities corresponding to the curve, a fitted stability curve in the plane is figured out, shaped like a parabola as in Fig. 5. In Fig. 5, the system stability curve divides the plane into two parts, the stable region which is shaded and unstable region which is white. Then for analyzing the relationship between MDP and quantitatively, the curve is fitted according to Connors' model since there are several empirical design guidelines based on the Connors' equation are developed [Hassan and Hossen (2010)]. In the engineering, the critical pitch flow velocity is always expressed as follows, where is the fluid density, and γ are empirical constants. is the structural natural frequency of tube in air circumstance and the simulations and referred experiment results in the paper are acquired in air flow. Over past five decades, a lot of experimental investigations are carried out in order to modify the constants to be suitable for different tube arrays. In this work, Connors' model is chosen as a reference to fit the stability curve and the pitch-diameter ratio chosen is 1.5. As shown in Fig. 5, the corresponding critical reduced flow velocities in 50 cases of MDP are computed. And several typical MDP and critical reduced flow velocities are listed in Tab. 2. Furthermore, the curve is fitted by these points in the form of Eq. (51) and values of and γ are obtained accordingly, which are 5.428 and 0.434 respectively. The fitted curve indicates that the critical flow velocity increases with MDP. In the lower MDP condition, it increases sharply and with the MDP increasing, the slope of the critical velocity growth becomes gentle continuously.   The Reynolds numbers depending on the mean gap flow velocity over the range of = 273 to = 1092 are listed in Tab. 3. The dimensionless flow velocity ( * = / ) and pressure distributions along the flow channel in a certain MDP under several flow velocities are also showed respectively in Fig. 6 and Fig. 7. As mentioned in the flow channel model, the locations along the fluid channel are described by one dimensional curvilinear coordinate , the velocity ( , ) and pressure ( , ) are computed along the coordinate , the distributions of which in transverse direction of flow channels are not presented. The location index ( * = / ) is introduced as shown in Fig. 6 and Fig. 7, which is the coordinate s scaled with the tube diameter. It respectively takes values of -1.5, 0, and 1.5 along the channel centerline at locations responding to the centers of tubes. It is observed that the velocity fluctuates along the coordinate , and it holds a peak level in the vicinity of the tube. And there is a general downward trend with the pressure distribution as it fluctuates along the coordinate in a certain Reynolds number. It shows that for low Reynolds numbers in scope of this paper where turbulence could be negligible, the velocity and pressure distribution respectively present similar trends along the channel coordinate while mean values rise with the increase of Reynolds number.  In this study, 50 cases of MDP respectively in 200 flow velocities are used to satisfy efficiency demand. However, discrepancies exist between the simulation curves and the scattered experimental data in Fig. 8. This indicates that due to the complexity of fluid-structural interaction, it is difficult to model and simulate the fluidelastic mechanism accurately by the models with much hypothesis and simplifications. Although the original flow channel model is progressive, it has deficiencies as a result of a series of assumptions and approximated parameters. The channel flow is assumed to be incompressible and inviscid, as well as wake phenomenon is ignored. Moreover, some artificial parameters such as the time lag, the position of separation and attachment points may also involve in errors. The fluid mechanics model is hence expected to be further improved in the future, which may approach the fluidelastic behavior more accurately.  Fig. 9. It demonstrates that the computation time is significantly reduced and efficiency is improved dramatically. In the Connors' Eq. (51), the tube instability only depends on MDP and the reduced flow velocity. The coupling effects between tube and fluid flow is not reflected. Now the impacts of fluid parameters will be discussed by the proposed method.
In the view of damping-controlled mechanism theory, the tube-flow coupling system includes both structural damping and fluid drag forces. The drag force coefficient ℎ related to the resistance influence during the process the fluid flow passes through the tube bundles, which can be calculated by solving the average pressure drop based on energy principle. In Fig. 10, the tube stability curve at different ℎ shows that, the greater stability region area appears with the value of ℎ increasing in − diagram. The reason is that higher ℎ value means more energy is dissipated for the coupling system, consequently the system has better stability in terms of energy conservation law. On the condition that the flow velocity is constant, in the SDM, the dimension of the transfer matrix is determined by the sample frequency, which also has an influence on the computation accuracy and efficiency of stability maps. The Fig. 11 shows the effect of sample frequency on the fitting curve coefficient and γ. It indicates that the parameters change significantly with the sample frequency when it is below 10 kHz. And when the sample frequency exceeds 10 kHz, the changes of parameters are small and negligible. Considering that the calculation amount also increases with sample frequency, 10 kHz is hence chosen as the sample frequency for the stability computation.

Summary
The proposed model is developed based on flow channel model and the tube motion is set to be arbitrary without limitation. The time delay between flow redistribution and tube vibration is also considered. The SDM is introduced in this study to analyze the stability problem of the dynamic equation with a constant delay term and a time delay integral term. The stability map of tube vibration is figured out visually based on eigenvalues of the built transfer matrix by SDM. Account for the unlimited tube motion and predicting the instability threshold directly via eigenvalues based on SDM, this proposed method presents definite improvement in computation efficiency with comparing to Hassan and Hayder's model. Numerical investigations of a single flexible tube within a rigid array subjected to crossflow are carried out for a typical square array. Effects of tube structural parameters and fluid related parameters may affect the stability are discussed. The stability depends on MDP mainly, which is determined by tube structural parameters. The stability is also influenced by fluid related parameters. When the fluid passes through tube bundles, greater fluid resistance may lead to more stable tube vibration than lower one. As the flow channel model is based on a series of assumptions and empirical parameters, in the future, a more elaborate fluid mechanics model is expected to be developed, which may predict the instability threshold with higher accuracy.