There exist wave motion mathematical models in infinitely long geometrically nonlinear shells filled with viscous incompressible liquid. They are based on related hydroelasticity problems, described by dynamics and viscous incompressible liquid equations in the form of generalized KdV equations. Mathematical models of wave process in infinitely long geometrically nonlinear coaxial cylindrical shells are obtained by means of the small parameter perturbation method. The problems differ from the already known ones by the consideration of viscous incompressible liquid presence. The system of generalized KdV equations is obtained on the basis of related hydroelasticity problems, described by shell dynamics and viscous incompressible liquid equations with corresponding boundary conditions. This paper deals with investigating of wave occurrence model of the two geometrically and physically nonlinear elastic coaxial cylindrical Kirchhoff – Love type shells containing viscous incompressible liquid between and inside them. The difference Crank – Nicholson type schemes aimed at investigating equations systems with the consideration of liquid impact are obtained with the help of Grobner basis construction. To generate these difference schemes, basic integral difference correlations, approximating the initial equations system, are used. The use of Grobner basis techniques makes it possible to generate the sch ¨ emes allowing to obtain discrete preservation laws analogues to the initial differential equations. To do this, equivalent transformations were made. On the basis of computational algorithm the software allowing to construct graphs and to obtain Cauchy problem numerical solution was developed, using the exact solutions of the coaxial shell dynamics equations system as an initial condition.