Based on the time-varying solution of layered fractional viscoelastic soils considering groundwater level, the interaction between soils and piles is studied by a finite element–boundary element coupling method. The single pile is regarded as a one-dimensional compressible bar to establish its stiffness matrix equation by the finite element method. Then, the fundamental solutions of layered fractional viscoelastic soils considering groundwater level are used as kernel functions in the boundary element method. According to the displacement continuity condition on the pile-soil interface, the interaction between single pile and soils is solved. Furthermore, by incorporating the pile-pile interaction as well as the balance and coordination conditions of the rigid cap, the solution of the single pile is extended to the pile group. Finally, the correctness of the presented method and program is validated, and several examples are designed to explore the effects of groundwater level, soil properties, pile parameters, and soil stratification on the pile-soil interaction. The analysis results show that, as the groundwater level declines, the displacement of the pile at the initial moment increases, and the reaction forces of center and side piles increase. Moreover, the increase of fractional order can accelerate the rate of pile displacement, reduce the axial force on the upper pile body, and decrease the reaction force at the top of the center and side piles.