We construct the exact solution of the quasi-static boundary value problem for a multilayer thick-walled tubes made of physically non-linear viscoelastic materials obeying the Rabotnov constitutive equation with two arbitrary material functions for each layer (a material creep compliance and a function which governs physical non-linearity). We suppose that every layer material is homogeneous, isotropic and incompressible and that a tube is loaded by time-dependent internal and external pressures (varying slowly enough to neglect inertia terms in the equilibrium equations) and that a plain strain state is realized, i.e. zero axial displacements are given on the edge cross sections of the tube. We obtained the closed form expressions for displacement, strain and stress fields via the single unknown function of time and integral operators involving this function, pairs of (arbitrary) material functions of each tube layer, preset pressure values and ratios of tube layers radii and derived integral equation to determine this unknown function. To derive it we split and solve the set of non-linear integral equations for unknown functions of time governing strain fields of every tube layer and for unknown interlayer normal stresses (depending on time). Assuming material functions are arbitrary, we proved that the total axial force at a tube cross section doesn’t depend on a number of layers, their thicknesses and pairs of material functions governing their mechanical behavior and on a history of loading although stresses and strains do. The axial force depends only on a tube radii and current values of given pressures. It proved to be equal to the axial force calculated for homogeneous linear elastic tube although axial stress depends on radial coordinate in the case of non-linear viscoelastic materials. Assuming the material functions that govern non-linearity of each layer material coincide with a power function with a positive exponent and assuming their relaxation moduli are proportional to a single (arbitrary) function of time, we constructed exact solution of the resolving functional equation, calculated all the integrals involved in the general representation for the tube stress field and reduced it to simple algebraic formulas convenient for analysis.