Convolution-based time-domain simulation for fluidelastic instability in tube arrays

A convolution-based numerical algorithm is presented for the time-domain analysis of fluidelastic instability in tube arrays, emphasizing in detail some key numerical issues involved in the time-domain simulation. The unit-step and unit-impulse response functions, as two elementary building blocks for the time-domain analysis, are interpreted systematically. An amplitude-dependent unit-step or unit-impulse response function is introduced to capture the main features of the nonlinear fluidelastic (FE) forces. Connections of these elementary functions with conventional frequency-domain unsteady FE force coefficients are discussed to facilitate the identification of model parameters. Due to the lack of a reliable method to directly identify the unit-step or unit-impulse response function, the response function is indirectly identified based on the unsteady FE force coefficients. However, the transient feature captured by the indirectly identified response function may not be consistent with the physical fluid-memory effects. A recursive function is derived for FE force simulation to reduce the computational cost of the convolution operation. Numerical examples of two tube arrays, containing both a single flexible tube and multiple flexible tubes, are provided to validate the fidelity of the time-domain simulation. It is proven that the present time-domain simulation can achieve the same level of accuracy as the frequency-domain simulation based on the unsteady FE force coefficients. The convolution-based time-domain simulation can be used to more accurately evaluate the integrity of tube arrays by considering various nonlinear effects and non-uniform flow conditions. However, the indirectly identified unit-step or unit-impulse response function may fail to capture the underlying discontinuity in the stability curve due to the prespecified expression for fluid-memory effects.

Abstract A convolution-based numerical algorithm is presented for the time-domain analysis of fluidelastic instability in tube arrays, emphasizing in detail some key numerical issues involved in the timedomain simulation. The unit-step and unit-impulse response functions, as two elementary building blocks for the time-domain analysis, are interpreted systematically. An amplitude-dependent unit-step or unitimpulse response function is introduced to capture the main features of the nonlinear fluidelastic (FE) forces. Connections of these elementary functions with conventional frequency-domain unsteady FE force coefficients are discussed to facilitate the identification of model parameters. Due to the lack of a reliable method to directly identify the unit-step or unit-impulse response function, the response function is indirectly identified based on the unsteady FE force coefficients. However, the transient feature captured by the indirectly identified response function may not be consistent with the physical fluid-memory effects. A recursive function is derived for FE force simulation to reduce the computational cost of the convolution operation. Numerical examples of two tube arrays, containing both a single flexible tube and multiple flexible tubes, are provided to validate the fidelity of the time-domain simulation. It is proven that the present time-domain simulation can achieve the same level of accuracy as the frequency-domain simulation based on the unsteady FE force coefficients. The convolution-based time-domain simulation can be used to more accurately evaluate the integrity of tube arrays by considering various nonlinear effects and non-uniform flow conditions. However, the indirectly identified unit-step or unit-impulse response function may fail to capture the underlying discontinuity in the stability curve due to the prespecified expression for fluid-memory effects.
Keywords Tube array Á Fluidelastic instability (FEI) Á Convolution Á Unit-step response function Á Unit-impulse response function Á Nonlinear fluidelastic (FE) force Nonlinear FE damping coefficient d(s) Dirac delta function q Fluid mass density n Tube damping ratio in flow n 0 Tube damping ratio in vacuum / Unit-step response function

Introduction
Long, slender tube arrays in heat exchangers are susceptible to various flow-induced vibrations, among which fluidelastic instability (FEI) has the greatest potential to induce short-term structural damage [1]. Due to the extremely destructive nature of this instability, considerable efforts have been made in this area in the last few decades [2][3][4][5][6][7][8]. It is now wellknown that there exist two distinct instability mechanisms, i.e. fluidelastic (FE) damping-controlled and FE stiffness-controlled mechanisms. Excellent reviews on this topic have been provided by several authors [9][10][11][12].
There is an increasing interest in advancing the time-domain analysis of FEI to more effectively account for any nonlinear effects (e.g. due to loosened supports, initial axial load, and FE force) and nonuniform flow conditions [13][14][15][16]. Both quasi-static and quasi-steady models have been widely used in time-domain analysis [17,18]. The basic assumption of a quasi-static model is that the FE force on a vibrating cylinder is equal to that on a stationary cylinder whose configuration is identical to the actual instantaneous configuration of the vibrating cylinder [19]. In contrast, a quasi-steady model assumes that the FE force on a vibrating cylinder is equal to that on a stationary cylinder whose configuration is identical to the actual instantaneous configuration modified by the instantaneous velocity [19][20][21]. Therefore, the FE forces of both models are determined by the instantaneous flow status, while any potential fluid-memory (unsteady) effects are neglected. The fluid-memory effects may significantly affect the energy transformation between a tube and its surrounding flow, and hence, both quasi-static and quasi-steady models are unable to predict FEI with satisfactory accuracy [12]. Although the quasi-steady model is often modified with a time delay, it is well-accepted that the time delay cannot fully capture the fluid-memory effects [22][23][24]. The unsteady model, which assumes that the FE forces on a vibrating tube are functions of its displacement, velocity, and acceleration, can satisfactorily capture the fluid-memory effects, while the model is expressed in a mixed time-frequency domain [25]. Detailed introduction and comparison of these theories can be found in, e.g. [12,19,26].
To effectively capture the fluid-memory effects in the time domain, the complex relationships between tube displacement/velocity and FE force can be simulated by a superposition of scaled and timeshifted fundamental responses, i.e. convolution. The unit-step response function and unit-impulse response function, as two fundamental functions in convolution-based time-domain analysis, have been widely used in the aerodynamics of plate-like structures, e.g. airfoils and bridge decks. Analytical solutions of the unit-step and unit-impulse response functions were made available for a thin airfoil [27]. However, for a bridge deck, a theoretical treatment of the unit-step or unit-impulse response function remains intractable due to complicated flow separation. As a result, the unit-step or unit-impulse response function of a bridge deck is often identified either experimentally through wind tunnel tests or numerically through computational fluid dynamics-based simulations. Both the direct method, in which the bridge deck is excited with a smoothed unit-step or unit-impulse displacement/velocity, and the indirect method, in which the function is numerically fitted according to the inverse Fourier transform of the frequencydependent aerodynamic transfer functions, have been advanced [28][29][30][31].
The application of a convolution-based time-domain simulation in flow-induced vibration analysis of tube arrays is very limited. Convolution was first introduced to simulate the FE forces of tube arrays by Granger and Païdoussis [22]. They modified the timedelayed term in a quasi-steady model with a unit-step response function-based convolution. The modified model (referred to as a quasi-unsteady model) showed a clearly improved predictive capability compared with the quasi-steady model. Using a vortex wake model, Meskell [23] showed that the unit-step response function can be approximated numerically without using experimental data. However, Granger and Païdoussis [22] and Meskell [23] both focused on the FEI of linear tube arrays in uniform flow, and hence, they converted the model into the frequency domain to determine the critical flow velocity. To the best of the authors' knowledge, the only convolutionbased time-domain analysis was conducted recently by Piteau [32], in which the linear FE force on a single flexible tube in a rigid tube array was simulated by a unit-step response function-based approach. Granger and Païdoussis [22] and Piteau et al. [32] identified the unit-step response function based on the vibration frequencies and damping ratios obtained in a free vibration test. Meskell [23] deduced the response function using a vortex wake model. These identification methods should be limited to a tube array with a single flexible tube. This paper presents a systematic interpretation of the unit-step response function and unit-impulse response function with a specific reference to the simulation of FE forces of tube arrays. Different from previous applications [22,23,32], the unit-step or unit-impulse response functions in this paper are determined based on the full unsteady FE force coefficients obtained from the forced vibration tests of Tanaka and Takahara [25]. It is assumed that the effects of each tube in an array on the FE forces of a tube can be summed linearly. Hence, the response functions determined based on FE force coefficients are applicable for an array with multiple flexible tubes. The connections of these elementary functions with frequency-domain unsteady FE force coefficients are introduced to facilitate the identification of model parameters. Some key numerical issues involved in the time-domain simulation, e.g. the integral of fluidmemory terms and the transient behaviour, are discussed. Amplitude-dependent unit-step and unitimpulse response functions are introduced to capture the main features of the nonlinear FE forces. Numerical examples of two tube arrays, containing both a single flexible tube and multiple flexible tubes, are provided to validate the fidelity of the time-domain simulation.
The subsequent parts of this paper are organized as follows: theoretical backgrounds of the unit-step and unit-impulse response functions, as well as their connections with frequency-domain unsteady FE force coefficients, are introduced in Sect. 2; the amplitude-dependent unit-step and unit-impulse response functions are introduced in Sect. 3; the indirect identification of unit-step and unit-impulse response functions are discussed in Sect. 4; a recursive function is deduced for FE force simulation in Sect. 5 to reduce the computational cost of the convolution operation; results of the numerical examples are provided in Sect. 6; the main conclusions are summarized in Sect. 7.
2 Unit-step and unit-impulse response functions for FE force

Unit-step response function
Considering an array of static tubes immersed in a two-dimensional flow, the flow-induced lift force acting on a tube can be expressed as where q is the fluid mass density; U is the pitch flow velocity; D is the diameter of the tube; and C L is the dimensionless lift coefficient of the static tube. If one of the tubes is disturbed by a unit-step motion in the cross-flow direction, the resulting lift force on the disturbed tube will exhibit a transient change described by where oC L =oY is the derivative of the static lift coefficient with respect to the dimensionless displacement Y = y/D in the cross-flow direction; y is the cross-flow tube displacement; / is the unit-step response function; and s is the dimensionless time. Inspired by the Jones approximation [33] of the analytical unit-step response function (i.e. the Wagner function [27]) of a thin airfoil, the following general expression is often employed in a case without an analytical solution: where a j and b j (j = 1 * n) are constants to be identified and b j [ 0; n is the order of the unit-step response function. Equations (2) and (3) indicate that in the case of the step motion, the lift force would, instead of entering a new steady state immediately, undergo a process of transient growth due to the fluid-memory effects.
The extension of Eq. (2) to a general tube motion, starting from rest at s = 0, can be obtained in terms of the well-known Duhamel's integral or equivalently, according to an integration by parts and a change of variables where prime represents the derivative with respect to dimensionless time s. Once the unit-step response function for a tube is available, the FE force induced by an arbitrary motion can be determined according to Eqs. (4) or (5).
Attempts have been made to directly measure the unit-step response function of a separated bluff body (e.g. a rectangular cylinder or a bridge deck) in wind tunnels or based on computational fluid dynamics simulations [28][29][30]. However, the direct identification method requires further improvement before it can become a routine identification technique. On the other hand, it is more practical to measure the unsteady FE force coefficients in the frequency domain, as pioneered by Tanaka et al. [25,34] and Chen et al. [26,35]. Therefore, a promising indirect method is to fit an equivalent unit-step response function based on the inverse Fourier transform of the measured unsteady FE force coefficients. To this end, Eq. (5) is transformed into the frequency domain as where the overbar represents the Fourier transform operator; K = 2pfD/U is the reduced frequency; and f is the tube vibration frequency in flow.
Noting that the Fourier transform of the exponential term e Àb j s in /(s) is 1/(b j ? iK), the bracketed terms in Eq. (6) can be expressed as On the other hand, the unsteady FE force on a tube is conventionally expressed in a mixed time-frequency domain as [25,34] FðtÞ or in a dimensionless form where c k and c d are the dimensionless FE stiffness and damping coefficients, respectively; the overdot represents the derivative with respect to time t. The equivalent expression of Eq. (9) in the frequency domain can be expressed as Comparing Eqs. (6) and (10), the relationships between the Fourier transform of the unit-step response function and the unsteady FE force coefficients may be given by It is then possible to identify the values of a j and b j based on c k and c d through least squares fit. Details of the identification procedure will be discussed later in Sect. 4.
For the in-line square tube array shown in Fig. 1, the FE forces on tube O are influenced by the vibrations of tube O and its surrounding tubes. To reduce the complexity of the FE force model, it is often assumed that only the effects of tubes R, L, U, and D need to be considered, and the effect of each tube can be summed linearly [25,26,34,35]. Following the variables defined by Tanaka et al. [25,34], the FE forces on tube O can be expressed as where j = O, R, L, D, and U denote the five tubes shown in Fig. 1; X is the dimensionless tube displacement in the streamwise direction; / YjY , / YjX , / XjX , and / XjY are unit-step response functions, in which the first, second, and third suffixes represent the direction of the force, the position of the vibrating tube, and the direction of the tube motion, respectively. As an example, / YRX represents the unit-step response function related to the y-direction force induced by the x-direction vibration of the right tube (R).

Unit-impulse response function
Considering again a static tube array immersed in a two-dimensional flow, if one of the tubes is disturbed by a unit-impulse motion in the cross-flow direction, the resulting lift will exhibit a transient change described by where I is the unit-impulse response function.
The FE force after a unit-impulse motion will exhibit a transient change from the force determined by the original steady position of the tube and finally fade out to the same value. Therefore, the unit-impulse response function can be divided into a steady part and an unsteady part where ½I s ð Þ s ¼ oC L oY d s ð Þ is the steady part, with d(s) representing the Dirac delta function; the unsteady part, ½I s ð Þ us , converges to zero with increasing time. The FE force on a tube due to an arbitrary motion can be expressed in terms of the unit-impulse response function as To identify the unit-impulse response function based on the unsteady FE force coefficients, Eq. (15) can be transformed into a frequency-domain expression given by The steady part can be expressed in the frequency domain as The frequency-domain expression of the unsteady part can be approximated by the well-known Roger's rational functions [36] Combining Eqs. (17) and (18), the frequencydomain expression of the unit-impulse response function is obtained as The equivalent time-domain expression of Eq. (19) is Comparing Eqs. (10) and (19), the relationships between the Fourier transform of the unit-impulse response function and the unsteady FE force coefficients are given by It is obvious from Eqs. (8), (11), and (21) that c j ¼ Àa j oC L oY and d j ¼ b j . Furthermore, Eqs. (5) and (15) indicate that the unit-impulse response function I(s) is the derivative of the unit-step response function / (s) with respect to dimensionless time s, with the exception of s = 0. Hence, the FE forces expressed in terms of unit-step and unit-impulse response functions are actually equivalent. This is exactly true as long as these functions are indirectly derived based on the frequency-domain unsteady force coefficients. However, when a direct identification method is utilized, a unit-impulse response function may be preferred over a unit-step response function [37]. If the surrounding tubes are positioned asymmetrically, the FE force due to a sudden upward unit-step motion may differ from that due to a sudden downward unit-step motion. In contrast, the unit-impulse response function should be able to capture the effects of this asymmetry since the tube is first disturbed in the upward direction and then downward to its original position.
It is well-known that the FE force can be wellapproximated by the quasi-steady theory if the reduced flow velocity U r = U/fD is sufficiently high (or equivalently, if the reduced frequency K is sufficiently low). According to the quasi-steady theory [38], the FE force can be given as where C D is the static drag coefficient. According to Eqs. (9) and (22), it is obvious that as K approaches zero, c k and c d should approach oC L oY and -C D , respectively. To ensure that the simulation of the unit-step or unit-impulse response function is identical to the quasi-steady theory as K approaches zero, it is necessary to constrain For the flexible tube array shown in Fig. 1, the FE forces on tube O can be expressed as FxðsÞ where I YjY , I YjX , I XjX , and I XjY are unit-step response functions. The suffixes of I are the same as those of the unit-step response function /.

Consideration of nonlinear FE force
The FE forces are essentially nonlinear functions of tube responses due to the complicated interactions between the fluid and the tubes. The nonlinear FE forces may result in limit cycle oscillations after the critical flow velocity for FEI [39] and/or lead to undesired large amplitude vibrations before the critical velocity (in a hysteresis region) [40]. Furthermore, the effects of nonlinear FE forces may be nonnegligible in evaluating the fretting wear or fatigue life of a tube array [41]. To achieve a more reliable structural dynamic analysis, there is a rising demand to consider nonlinear FE forces. Some attempts have been made to consider the nonlinear FE effects utilizing the quasisteady model [39] and the semi-analytical flow channel model [40,41]. However, these models should be used with caution since they include assumptions that might be unrealistic.
The unsteady model can also be modified with higher-order nonlinear damping and/or stiffness terms to consider the nonlinear FE forces in the frequency domain [42,43]. As an example, Meskell and Eret [42] introduced a cubic damping term to simulate the postcritical limit cycle oscillations of a flexible tube surrounded by a rigid normal triangular tube array: where b is a nonlinear FE force coefficient.
To consider the nonlinear effects in the time domain, the Volterra theory [44] provides a promising choice as a direct extension of the linear convolution. However, the identification of higher-order kernels in a Volterra-based model is challenging, and the computational cost increases exponentially with increasing model order. On the other hand, it is known that both the postcritical limit cycle oscillation and hysteresis behaviour of a FE system are essentially related to the amplitude dependency of FE forces [45,46]. Inada et al. [47] highlighted the amplitude dependency of frequency-domain FE force coefficients for an in-line square tube array and further proposed predicting postcritical vibrations with amplitude-dependent FE force coefficients. Accordingly, it might be beneficial to develop an amplitude-dependent unit-step or unit-impulse response function to capture the main features of the nonlinear FE forces in the time domain. For example, the relationship between the amplitude-dependent unit-step response function /(A, s) and frequency-domain FE force coefficients are given by where A is the vibration amplitude. The amplitudedependent FE force coefficients, c k (K, A) and c d (K, A), can be identified based either on the displacement signals of free vibration experiments [42], or on harmonic forced vibration experiments with various vibration amplitudes [43]. For a flexible tube array, the nonlinear FE forces on a tube, e.g. tube O in Fig. 1, are then approximated as FxðsÞ A practical implementation of Eq. (26) is shown in Sect. 6. The FE force can also be expressed in terms of amplitude-dependent unit-impulse response functions, which is not repeated for brevity. It is noted that the physical meanings of the nonlinear amplitude-dependent unit step/impulse response functions are uncertain. These functions are merely time-domain equivalents of the amplitude-dependent, unsteady FE force coefficients in the considered range of reduced velocities. Despite the lack of a sound physical meaning, it will be shown that the amplitude-dependent unit step/impulse response function is capable of capturing the main nonlinear feature (i.e. amplitudedependency) of the unsteady FE forces. Since the amplitude-dependent unit step/impulse response function is identified based on the unsteady FE force coefficients (obtained from forced vibration data), the present nonlinear model is believed to be more accurate than the quasi-steady model in calculating the large-amplitude vibrations in the post-instability regime.

Identification of unit-step and unit-impulse response functions
The values of a j and b j (or c j and d j ) can be obtained through nonlinear least squares fit based on Eq. (11) (or Eq. (21)). As an example, according to Eq. (11), the values of a j and b j can be obtained by minimizing the following error function where Z represents the sampling number of the unsteady FE force coefficients involved in the identification; w 1 and w 2 are the weight coefficients of the FE stiffness and damping parts, respectively.
On the other hand, it might be more reasonable to calculate these values to minimize the fitting errors between the forces simulated by the unsteady force coefficients (Eq. (9)) and the forces simulated by the convolutions (Eqs. (4) or (15)). As noted by Piteau [32], the apparently large errors of the fitted c d and c k at a low reduced flow velocity do not introduce significant errors in the calculated FE responses. Indeed, according to Eq. (9), the vibration frequency and damping ratio of the tube at any reduced frequency can be deduced as where m is the tube mass plus the fluid added mass per unit tube length; f 0 is the tube natural frequency in a vacuum; n 0 is the tube damping ratio in a vacuum; and n is the tube damping ratio in a flow. It is then obvious that the FE stiffness and FE damping forces are proportional to c k /K 2 and c d /K, respectively. Hence, the parameters of the unit-step response function can be obtained by minimizing the following error function Similarly, the parameters of the unit-impulse response function can be obtained by minimizing As presented in Eq. (29a) (or Eq. (29b)), both a j (and c j ), located in the numerators, and b j (or d j ), located in the denominators, are unknown parameters to be identified. Conventionally, this nonlinear fitting problem involves selecting a group of initial values and an iterative numerical procedure to find a local optimal solution. The results of such methods are sensitive to the initial values, and sometimes it might be necessary to repeat the whole process by starting with different groups of initial values.
To this end, the nonlinear fitting problem can be degenerated to a linear fitting problem by searching b j (or d j ) within a predefined range. The identification procedure of the unit-step response function is illustrated in Fig. 2 as an example. For specific values of b j , the optimal values of a j can be conveniently identified through a linear fitting. Hence, the nonlinear fitting is linearized during each searching step, and the procedure leads to the optimal solution in the searching range (b min , b max ) rather than a local optimal solution that is often obtained in a nonlinear fitting. The unitimpulse response function can be identified with a similar procedure, which is not repeated for brevity.
It is important to note that the unit-step or unitimpulse response function identified based on the unsteady FE force coefficients may not physically represent the true properties of the fluid-structure interaction. Hence, the function is not expected to be applicable outside the range of reduced velocities corresponding to the original experimental data.

Integral of fluid-memory terms
A time-domain analysis of FEI entails a step-by-step process in which Eqs. (4) or (15) is employed to simulate the FE force. In both expressions, there are fluid-memory terms related to the whole time-history of the tube motion, which may be extremely timeconsuming in terms of the convolution integration. The significant computational cost might prevent the practical application of the time-domain simulation in cases of three-dimensional FE analyses of flexible tube arrays. Fortunately, the computational cost can be largely reduced if the fluid-memory terms are transformed into recursive functions. Taking the unit-step response function as an example, Eq. (4) can be reformulated as where The first term on the right-hand side of the last equal sign in Eq. (30), i.e. 1 2 qU 2 D oC L oY Y s ð Þ, is only dependent on the displacement of the current instance.
However, the second term, i.e. P n j¼1 A j s ð Þ, depends on the entire time history of the tube motion. A recursive function can be generated if there exists an explicit relationship between A j at dimensionless time step s ? Ds, and its value at dimensionless time step s. A j at dimensionless time step s ? Ds can be expressed as Assuming the integrand of the second term on the right-hand side of the last equal sign in Eq. (32) changes linearly between each dimensionless time step, it might be obtained as It is then possible to calculate the FE force F(s ? Ds) recursively based on F(s) according to Eqs. (30) and (33). A similar recursive function can be deduced for the unit-impulse response function-based convolution, which is not repeated for brevity. It should be highlighted that to consider the amplitude dependency of FE forces, a j and b j should be functions of the vibration amplitude, as presented in Eq. (25). Hence, in a time-domain simulation considering nonlinear FE effects, it is necessary to update a j and b j at each step according to the instantaneous vibration amplitude.

A linear case
A numerical example of an in-line square tube array with pitch ratio P/D = 1.33 is provided to verify the fidelity of convolution-based time-domain simulation. The frequency-domain FE force coefficients were measured by Tanaka et al. [25], as given in Fig. 3. The unit-step and unit-impulse response functions are fitted based on these coefficients. Second-order functions (i.e. n = 2) are adopted for / yOy , / yLy , / yLx , and / xLy (or I yOy , I yLy , I yLx , and I xLy ), while n = 3 is adopted for / xOx and / xLx (or I xOx and I xLx ) to increase the fitting accuracy. Further increasing the model order does not significantly increase the fitting accuracy. The limiting values of the searching range are b min = 0.01 and b max = 4.0. Parameters a j and b j of the resulting unit-step response functions are listed in Table 1. Parameters c j and d j of the resulting unitimpulse response functions are c j ¼ Àa j It is noted that both unit-step and unit-impulse response functions show slight overshooting features. The overshooting feature has also been observed for some separated bluff bodies (e.g. a rectangular cylinder or a bridge deck) [28,30] and for a flexible tube in a rigid triangular tube array by Li and Mureithi [24]. Generally, the overshooting phenomenon is believed to be induced by (or at least accompanied by) the flow separation from the structure. Due to the complicated interaction between fluid and structure, it is possible that the transient FE force exceeds the steady-stead force, which is reflected as the overshooting phenomenon in the response function. However, it will be shown later that the overshooting of an indirectly identified response function is quite arbitrary and might simply be a result of numerical fitting.
/ yOy and / xLy (or I yOy and I xLy ) approach to the unit (or zero) after 70 and 200 dimensionless time steps, respectively. The fading time indicates that after a unit disturbance of tube O or tube L, the resulting FE forces on tube O do not become stable until the flow transports 70D or 200D downstream. The tremendously long fading time is absolutely not physical and is a result of the numerical fitting. The fading time can be reduced by searching b j in a (b min , b max ) range with a larger b min . Indeed, if a function with shorter fading time is used, the fitted c d and c k will approach their quasi-steady values (hence, c d Á U r or c k Á U 2 r versus U r will approach a linear or quadratic curve) at a lower reduced wind speed U r . However, most of the experimental data (especially c d, yLx ) in Fig. 3 do not show such a trend, which might be due to a systematic experimental error, as discussed by Païdoussis et al. (2010).
Once the unit-step or unit-impulse response function is available, the frequency-domain FE force coefficients can be calculated based on Eqs. (11) or (21) to verify the fidelity of the indirect identification procedure discussed in Sect. 4. The recalculated FE force coefficients are presented in Fig. 3 together with the experimental measurements. The comparisons demonstrate that numerical approximations provide satisfactory results. Further increasing the function order n can produce slightly improved fittings, while the computational costs will also be increased. Hence, unit-step or unit-impulse response functions with n = 2 or 3 will be employed in the subsequent timedomain simulations.
The critical pitch flow velocity U cr can be calculated in the frequency domain according to the FE force coefficients (through a complex eigenvalue analysis) or in the time domain according to the unit-step or unit-impulse response functions. To damping ratio is n 0 = 0.2%, and the mass-damping parameter 2pmn 0 /qD 2 is varied by changing the mass ratio m/qD 2 . The resulting U cr /f 0 D are plotted versus the mass-damping parameter 2pmn 0 /qD 2 in Fig. 6. The stability boundaries obtained by the unit-step and unit-response functions are identical, and hence, only one curve is shown for the time-domain simulations of each configuration. The comparisons prove that the time-domain simulations can satisfactorily calculate the critical state of FEI. It is stated that the nonlinearities (from structural and FE aspects) and nonuniform flow conditions can be implemented straightforwardly in the time-domain simulations. On the other hand, to consider these nonlinear and/or non-uniform effects in the frequency domain, an iterative process is often required, and the computational cost of iteration increases significantly as the nonlinear and/or nonuniform effects become strong [19,32]. As noted in Fig. 6, there is a discontinuity in the stability curve of the frequency-domain simulation between 1 \ 2pmn 0 /qD 2 \ 10. A discontinuity in the stability curve of a tube array has also been noticed by Chen [3] and Price [12]. Hence, it is believed that the discontinuity is an actual feature of the FE system. However, this discontinuity does not exist in the curve of the time-domain simulation. The inconsistency is a result of numerical fitting in the indirect identification of unit-step or unit-impulse response functions. Due to the prescribed forms of these response functions, the variations of the captured FE stiffness and damping effects versus reduced flow velocity follows the rule specified by the response functions. As a result, the captured FE stiffness and damping effects (and hence the predicted stability boundaries) are actually smoothed compared with those of the original FE force coefficients. This is a shortcoming of the timedomain simulation, which is difficult to be improved by simply increasing the order of the unit-step or unitimpulse response function. This is because the frequency-domain equivalent of a unit-step or unitimpulse function should be a continuous function of the reduced frequency (see Eqs. (11) or (21)) regardless of the function order n, which makes the response function unable to capture the discontinuity in the FE effects. For the considered tube configuration, only one stability region exists. Multiple stability regions may exist for other configurations, e.g. a tube row [48], while the capability of the convolution-based time-   The fidelity of the time-domain simulation is further verified by comparing the simulated vibration frequency and damping ratio in the pre-instability and post-instability ranges. This is done for the first configuration (i.e. a single flexible tube in a rigid tube array) with a mass-damping parameter 2pmn 0 /qD 2-= 20, as presented in Fig. 7. It is noted that the assumption of linear superposition no longer holds in the post-instability range due to the significantly increased vibration amplitude. The vibration frequency and damping ratio at a specific flow velocity can vary remarkably with increasing vibration amplitude [49], which cannot be captured by a linear FE force model. The comparison in Fig. 7 proves that the time-domain simulation can predict not only the critical flow velocity but also the pre-instability behaviour of the tube. The vibration frequency increases with increasing U r due to the positive FE stiffness effects (reflected by the negative c k, yOy values in Fig. 3b), while the damping ratio decreases due to the negative FE damping effects (reflected by the positive c d, yOy values in Fig. 3a). Furthermore, the simulation can be extended into a nonlinear convolution scheme to consider the nonlinear FE force in the post-instability range, which will be discussed in the next example.

A nonlinear case
A numerical example of a single flexible tube in a normal triangular tube array with pitch ratio P/ D = 1.33 is provided to prove the fidelity of timedomain simulation based on amplitude-dependent unit-step or unit-impulse response functions. Postcritical limit cycle oscillations of this configuration were experimentally measured by Meskell and Fitzpatrick [16]. Meskell and Eret [42] further reported the FE force coefficients (i.e. c k , c d , and b, as shown in Eq. (24)) identified from the experimental displacement responses. The amplitude-dependent FE force coefficients c d, yOy can be reconstructed from c d and b according to an averaging method [50,51], and the results at vibration amplitudes A = 0, 0.002, and 0.004 mm are given in Fig. 8a. Figure 8b presents the results of c k, yOy , which are almost independent of vibration amplitude according to [16,42].
Second-order unit-step response functions are fitted based on these coefficients. Parameters a j and b j of the resulting unit-step response functions are listed in Table 2. FE force coefficients fitted by the amplitudedependent unit-step response functions are presented in Fig. 8 together with the experimental measurements. The comparisons demonstrate that amplitude-dependent unit-step response functions can satisfactorily capture the amplitude dependency of FE force coefficients.
The postcritical vibration amplitude can be calculated in the frequency domain according to the amplitude-dependent FE force coefficients (through an iterative process) or in the time domain according to the amplitude-dependent unit-step response functions. To compare the results of time-domain and frequency-domain simulations, the postcritical vibration amplitudes are calculated for this specific configuration with the following mechanical properties: n 0 = 1.5%, f 0 = 6.5 Hz, and 2pmn 0 /qD 2 = 177.4. In the frequency-domain simulations, the vibration amplitudes are calculated by a regular iterative procedure . The resulting postcritical vibration amplitudes are plotted versus pitch flow velocity in Fig. 9. The comparisons prove that the time-domain simulation can satisfactorily calculate the postcritical response.

Transient behaviour of time-domain simulation
Since the unit-step or unit-impulse response function is identified indirectly based on the frequency-domain FE force coefficients, it should be remembered that they carry no additional information than that contained in these FE force coefficients. As a result, the shortcomings of conventional frequency-domain analysis, such as the limited applicable reduced frequency range, are remained in the time-domain analysis. In  Fig. 9 Post-critical responses of a flexible tube in a rigid normal triangular tube array addition, extra uncertainties can be introduced by numerical fitting. Hence, the indirectly identified unitstep or unit-impulse response function may not correctly reflect the physical fluid-memory effects in terms of both magnitude and duration. It is natural to expect that the distortion of true fluid-memory effects can be reduced if FE force coefficients in a wider reduced frequency range are employed in the identification. However, these coefficients are often only available in a limited reduced frequency range, especially if they are identified in a free vibration test based on the variations of the tube vibration frequency and damping ratio versus flow velocity [32,52]. In such a case, significant distortion of the fluid-memory effects can be introduced by inappropriate choices of b j or d j , for the unit-step or unit-impulse response function, respectively.
To demonstrate the possible distortion of the fluidmemory effects captured by an indirectly identified function, FE force coefficients c d, yOy and c k, yOy of the in-line square tube array within the reduced flow velocity range 3.5 B U r B 25 are utilized to identify the unit-step response function. Three independent runs of the identification procedure were carried out, in which b j values were searched within (0.08, 4), (0.001, 4), and (0, 4), respectively. The resulting unitstep response functions, denoted as / 1 , / 2 , and / 3 , respectively, are given in Fig. 10. Table 3 lists the a j and b j values of the three unitstep response functions. The FE force coefficients fitted by three unit-step response functions are shown in Fig. 11 together with the original experimental measurements. The unit-step response functions all produce reasonable fittings of the FE force coefficients. However, the identified unit-step function is largely dependent on the selection of the search range.   the unit-step response function, as shown by / 2 and / 3 in Fig. 10. The FE force induced by a harmonic motion for U r = 20 is simulated by the unit-step response functions, as presented in Fig. 12, together with that simulated by the frequency-domain FE force coefficients. It is noted that, due to the long fading time of the unit-step response functions, the FE forces simulated by / 2 and / 3 depart remarkably from those simulated by FE force coefficients. After a short evolution, the force simulated by / 2 collapses to the frequency-domain result. The force simulated by / 3 will eventually match the frequency-domain result after sufficient evolution. However, this cannot be seen in Fig. 12 since the evolution time is very long (approximately 6500 dimensionless time steps). Therefore, it is necessary to assign a reasonable range for b j or d j in the indirect identification of the unit-step or unit-impulse response function.
Since all three unit-step response functions in Fig. 10 are fitted based on the unsteady FE force coefficients, it is not possible to tell which of the three functions is closest to the true physical problem. Indeed, none of the three response functions are believed to be physically true. To obtain a response function that reflects the true physical problem, transient information (e.g. forced vibration data with unit-step input or amplitude-varying input) is required for the identification, which is not within the scope of this paper. In addition, the values of b min and b max are determined empirically in this paper. Fortunately, the accuracy of the response function (in terms of simulating the FE force) is insignificantly affected by the searching range as long as b min is not too small.

Conclusions
This paper presents a detailed interpretation of a convolution-based numerical algorithm for the timedomain analysis of the fluidelastic instability (FEI) of tube arrays. The unit-step and unit-impulse response functions, as two elementary building blocks for the convolution-based time-domain simulation, are interpreted systematically. An amplitude-dependent unitstep or unit-impulse response function is introduced to capture the main features of the nonlinear fluidelastic (FE) forces in the time domain. Theoretically, there is a close relationship between these elementary functions and the conventional frequency-domain unsteady FE force coefficients. Due to the lack of a reliable method to directly identify the unit-step or unit-impulse response function, it is more realistic to identify the function indirectly based on the FE force coefficients. The indirectly identified unit-step or unitimpulse response function carries no additional information than that contained in the FE force coefficients, which are often only available in a limited reduced frequency range. Hence, the indirectly identified function is not physically true. However, both functions may play an effective role in time-domain FE analysis since they correctly capture the frequency spectrum of FE forces.
Numerical examples of two tube arrays, containing both a single flexible tube and multiple flexible tubes, are provided to validate the fidelity of the time-domain simulation. Numerical results prove that the timedomain simulation can achieve the same level of accuracy as the frequency-domain simulation based on the unsteady FE force coefficients.. Time-domain analysis has the superiority of conveniently implementing any nonlinearities (originating from structural and FE aspects) and nonuniform flow conditions. However, the indirectly identified unit-step or unitimpulse response function can fail to capture the discontinuity in the stability curve since the variations of FE stiffness and damping versus reduced flow velocity are smoothed by numerical fitting.
Funding Open access funding provided by NTNU Norwegian University of Science and Technology (incl St. Olavs Hospital -Trondheim University Hospital). This study was funded by the Research Council of Norway (strategic Norwegian research program PETROMAKS2, grant agreement number 280713).

Declaration
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.