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.