The aim of the paper is to obtain a system of nonlinear evolution equations for two coaxial cylindrical shells containing viscous fluid between them and in the inner shell, as well as numerical modeling of the propagation processes for nonlinear solitary longitudinal strain waves in these shells. The case when the stress-strain coupling law for the shell material has a hardening combined nonlinearity in the form of a function with fractional exponent and a quadratic function is considered. Methods. To formulate the problem of shell hydroelasticity, the Lagrangian–Eulerian approach for recording the equations of dynamics and boundary conditions is used. The multiscale perturbation method is applied to analyze the formulated problem. As a result of asymptotic analysis, a system of two evolution equations, which are generalized Schamel– Korteweg– de Vries equations, is obtained, and it is shown that, in general, the system requires numerical investigation. The new difference scheme obtained using the Grobner basis technique is proposed to discretize the system of evolution equations. Results. The exact solution of the system of evolution equations for the special case of no fluid in the inner shell is found. Numerical modeling has shown that in the absence of fluid in the inner shell, the solitary deformation waves have supersonic velocity. In addition, for the above case, it was found that the strain waves in the shells retain their velocity and amplitude after interaction, i.e., they are solitons. On the other hand, calculations have shown that in the presence of a viscous fluid in the inner shell, attenuation of strain solitons is observed, and their propagation velocity becomes subsonic.