On the mathematical models of the Timoshenko-type multi-layer flexible orthotropic shells

Mathematical models of multi-layer orthotropic shells were reconsidered based on the Timoshenko hypothesis. A new mathematical model with ε\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon $$\end{document}-regularisation was proposed, and the theorem regarding the existence of a generalised solution to the model was formulated and proved. The algorithms of numerical investigation of models studied with the aid of the variational-difference method were developed. The associated stability problem was also addressed. A comparison of the results yielded by the considered models was carried out and discussed for numerous factors and parameters.


Introduction
The first papers devoted to the investigation of vibrations of piecewise non-homogeneous (with respect to thickness) multi-layer/composite plates and shells with the account of the geometric nonlinearity appeared in the sixties of the twentieth century. In these works, numerical computations of three-layer structural members with lightly stiffened internal layer (filler) have been used to show that inclusion of transverse shear and normal forces as well as deformations in the filler is required during modelling. Usually, owing to the complexity of nonlinear PDEs governing the dynamics of multi-layer plates and shells made from either orthotropic or anisotropic composite materials, the solutions to PDEs have been constructed by taking into account only small nonlinear vibrations and the firstorder approximation.
The first pioneering results related to the investigation of three-layer plates and shells have been reported by [1][2][3][4][5][6][7][8], and other. In references [1,2], three-layer orthotropic rectangular plates subjected to normal loads have been studied. Within the study, the governing PDEs have been derived and the dependency of nonlinear vibrations on the amplitude and frequency has been obtained in the first approximation.
Yu [5] found the approximate solution to the problem of nonlinear axially symmetric vibrations of a three-layer closed circular cylindrical shell with a light non-compressed filler and the membrane external layers. Kholod [7,8] has utilised the same deformation hypotheses regarding the filler while solving the problem of nonlinear transverse vibrations of a three-layer cylindrical simply supported panel symmetric with respect to its height, having a stiff filler and carrying loads by external layers.
The influence of deformations of the transverse shear on nonlinear vibrations of composite plates and shells has been studied in a vast number of references such as [6,[9][10][11][12]. In one of these works [6], it was pointed out that the influence of the transverse shear cannot be omitted while investigating nonlinear vibrations and stability of freely supported multi-layer plates, and it plays a crucial role in studying clamped multi-layer plates.
Kogan and Yurchenko [12] studied the influence of geometrical and mechanical parameters (relative plate thickness h/b, slenderness l = a/b, stiffness of the filler against the shear factor, bending stiffness of the carrying load layers, parameter e responsible for external damping of vibrations, the number of half-waves m and n) as well as of the type of boundary conditions on the character of free/forced nonlinear vibrations and the shape of the amplitude-frequency bone curves for a three-layer rectangular plate with a nonsymmetric (with respect to thickness) structure with a stiff, transversally isotropic filler transmitting the transverse shear and two isotropic load-carrying layers.
In general, factors such as a multi-layer mechanical structure, the presence of layers with different physical and geometrical characteristics, interlayer contact conditions, well-posed statements of the multi-layer strain as well as proper implementation of boundary conditions require more advanced modelling due to phenomena that are more complex than those exhibited by simple, reduced-order models.
On the other hand, mathematical models derived in the form of coupled nonlinear PDEs require the employment of more rigorous mathematical tools to reduce infinite dimensional problems to finite dimen-sional ones. In other words, it is necessary to reduce PDEs to truncated systems of nonlinear ODEs, putting emphasis on adequate modelling of the rich nonlinear structural behaviour.
Particular features of numerous directions of the theory development and the construction of mathematical models of multi-layer plates and shells have been described in references [13,14]. In general, one can distinguish two main approaches to investigating nonlinear vibrations of multi-layer plates and shells. Namely, one can employ: (i) a hypothesis of a straight, nondeformed normal (to all layers), and reduce the problem to computing a quasi-homogeneous plate/shell with reduced elastic parameters [15][16][17][18][19][20][21][22]; (ii) non-classical 2D theories based on integral hypotheses for the entire packet of layers, taking into account deformations and stresses of the transverse shear in all layers [5,6,[23][24][25][26][27][28].
Ugrimov [29] presented the layer-wise generalised theory of multi-layer shells based on the expansion of the displacement vector components of each layer into power series of the transverse coordinate. The results were validated by the analysis of the strain-stress states of one-and three-layer structures.
Fernandes [30] proposed a so-called mixed formulation to study elastic multi-layer plates for which the continuity of transverse shear stress is satisfied in a natural way. The introduced model was validated by bilayer and sandwich plates.
Basar and Eckstein [31] developed a finite multilayer shell element formulation with the account of large inelastic strains and finite rotations. The derived model was enforced into a multi-layer shell kinematics, putting emphasis on the proper modelling of stresses across the thickness. Also, a few examples validating the introduced theory and computational algorithms were reported.
Recently, Amabili and his co-workers have constructed a series of geometric nonlinear theories of the third order for the shells of general forms. For this purpose, they have employed five, seven, or eight parameters [34][35][36]. Since the terms of a higher order with respect to the transverse coordinate are also taken into account, the theory is feasible also for thick shells. The developed theory allows for investigating largeamplitude vibration excited by periodic load and the obtained results coincide with those yielded by other theories. Owing to complexity of the analysed problems, the investigations have been carried out with different numerical methods, i.e. the finite difference method, the finite element method, and the Galerkin methods.
Gutiérrez Rivera et al. [37] developed a new 12parameter finite element method of the shell in order to analyse large deformations of composite shell structures by using the third-order model.
In reference [38], the finite element method has been successfully used to study laminated composite skewed hyper shell roof, taking the account of an oblique impact with friction boundary elements.
Tornabene et al. [39] studied a mechanical behaviour of damaged laminated composite plates and shells by employing the higher-order shear deformation theories and using the accurate and stable generalised differential quadrature method (GDQ).
Large vibrations of three-and multi-layer plates, taking into account the initial imperfection, have been studied in references [12,[24][25][26]. In particular, in references [24,25] the parabolic form of distribution of deformation of the transverse shear along the thickness of the whole layer has been investigated. Five input differential equations have been reduced to one nonlinear second-order equation with square and cubic nonlinearity. Then, the equation has been solved using the multiscale approach. It has been also proved that in the case of either strong or weak nonlinearity, the behaviour of the plate essentially depends on the initially introduced imperfections [25].
More recently, Wang et al. [32] studied spherical cloaking by using multi-layer ordinary dielectric materials and optimised the total scattering cross section of a spherical multi-layer shell with a metallic core.
Jiang et al. [33] implemented a reduced-order model to evaluate the dynamics response of multi-layer plates under impulsive loads. The introduced approach aimed at either absorbing or transmitting the energy in transverse directions is based on the reverberation matrix method matched with the theory of generalised rays. The dynamic response of multi-layer plates was quantified using the generalised ray theory and the inverse Fourier transformation.
The authors of the present paper have been investigating the problems of statics and dynamics of multilayer structural members for many years. In particular, the multi-layer flexible Euler-Bernoulli and Timoshenko beams were investigated in reference [40], whereas multi-layer beams with a clearance have been studied in references [41,42]. Also, the non-classical mathematical models of coupled problems of thermoelasticity were derived and studied for multi-layer shells with initial imperfections in [43].
This paper extends the former investigations. It is organised in the following way. Section 2 comprises the description of fundamental hypotheses and relations as well as the derivation of the governing equations. Section 3 presents the modification of the Timoshenko model and the ε-regularisation model. It includes rigorous mathematical considerations including formulation of the Theorem 1 and its proof. Numerical investigations of stability of the multi-layer orthotropic shells within the earlier introduced improved theories are carried out in Sect. 4. Section 5 reports a comparison of the stability curves of the multi-layer symmetric shells based on the "load deflection" relations yielded by the four studied different shell mathematical models. Finally, Sect. 6 summarises the carried-out research.

Fundamental hypotheses
A multi-layer shallow shell being a part of the 3D space R 3 is considered. The rectangular coordinate system of the space R 3 is introduced as follows. An arbitrary surface of the shell, called a reference surface and defined as x 3 = 0, is chosen; the axes 0x 1 ,0x 2 are directed along the main curvatures of this surface, whereas the axis 0x 3 is directed towards the reference surface curvature axis. In the given system of coordinates, the shell is defined in the following manner: stands for the rectangular shell planform; δ n+m − δ 0 = 2h 0constant shell thickness; x 3 = δ 0 − -lower shell surface; x 3 = δ n+m − -upper shell surface;shell thickness measured from the upper surface to the reference surface x 3 = 0.
To develop mathematical models (MM2) based on the Timoshenko hypothesis, the following assumptions regarding the shell construction, material properties, layers and exploitation conditions are introduced.
1. The thickness of the i-th shell layer is denoted by δ i , i = 0, . . . , n + m the; m stands for the number of layers up to the layer consisting the surface x 3 = 0; n is the number of the remaining layers.
The interval x 3 ∈ (δ 0 − , δ n+m − ) is divided into subintervals with respect to x 3 within one layer 2. It is assumed that the shell is fully shallow, i.e. the coefficients of the first squared form of the reference surface are taken as A i ≈ 1. The main curvatures and the twist curvature of x 3 = 0 are constant (the Gauss curvature is neglected). Furthermore, it is assumed that the total thickness of layers is small compared to the curvature radii, i.e. 1 + K i j x 3 ≈ 1. 3. It is assumed that normal stresses σ 33 = 0 are small compared to the stresses in the state equations. 4. It is assumed that the shell is under a stationary temperature field. The temperature distribution is yielded by a solution to the Fourier heat transfer equation. If T 0 stands for the input shell temperature and T (x 1 , defines the grad as a vector of heat flow. 5. The considered piecewise homogeneous shell is composed of an arbitrary number of layers of the same thickness, but different stiffness. The layers are arbitrarily located with respect to the reference surface x 3 = 0. If the layers are orthotropic, then in each of them, there is one plane of elastic symmetry which is parallel to the plane tangent to the reference surface, and the two remaining planes are perpendicular to the axes Ox 1 , Ox 2 . 6. The deformable shell state is considered based on the assumption that deflections of the points of the reference surface can be of the same order as the shell thickness [44]. 7. The components of the displacement vector of the surfaces x 3 = 0 are denoted by are displacements of an arbitrary shell point, and u z 3 = u 3 (x 1 , x 2 , x 3 ) is the shell deflection at an arbitrary point. 8. The external shell surfaces are subjected to the load vector of the form: q + = q;q = 0, i.e. the shell behaviour caused only by the action of the transverse force is studied. 9. Since only thin flexible shells are considered, following monograph [9], the deformations in an arbitrary shell layer have the following form 10. One can transit from 3D to 2D theory of shells by taking into account the hypothesis of the second approximation, i.e. the Timoshenko-type hypothesis ( Fig. 1).
It is assumed that tangential displacements u z 1 , u z 2 , u z 3 are distributed along the shell thickness in a way described by the following linear form denote angles of the normal rotation with respect to the surface x 3 = 0 generated by deformations in the planes x 2 0x 3 and x 2 0x 3 , respectively. Then, equation (1) with the account of equation (2) yields Now, if one denotes by ε 11 , ε 22 , ε 12 the tangential deformations of the middle surface as follows the bending deformations, and by The Duhamel-Neuman principle [45] employed to each orthotropic (i-th) layer results in the following formulas The first three equations of (6), taking into account (8), yield where: E i 1 , E i 2 , E i 3 are the elasticity moduli; ν i 12 , ν i 13 , ν i 32 , ν i 31 , ν i 23 -Poisson's coefficients; G i 12 , G i 13 , G i 23shear moduli of the i-th layer. For the chosen form of orthotropy of the material layers, the following relations hold [5]: Now, one can take into account the static hypothesis σ i 33 = 0, find a formula for deformations ε i 33 , and substitute it into the first two relations of (6) using (8), (11). Then, in each shell layer, stress tensor components take the following form: Introducing the notation and taking into account (4), one can recast the stress formula to the following form suitable for further considerations: The equilibrium equations and the boundary conditions are obtained using the principle of virtual displacements [47]. This implies δV = δ A, where V stands for the potential energy of the shell deformation, and A is the work done by the external forces. One gets To obtain 2D equations of the shell equilibrium, the series (4) and the rule for computing of the triple integrals are used. For each layer, the interval with respect to Then, the notation n can be used and, based on the analogy to the one-layer shell, one denotes by T 11 , T 22 , T 12 the internal stresses; by Q 1 , Q 2 -transverse forces; byM 11 , M 22bending moments, and byM 12 -torsional moment: and hence The stress function stands for F (x 1 , x 2 ) and (14) is substituted into (15), where: The coefficients of the series (17) take the following form: and, owing to (10), It means that C 12 = C 21 , K 12 = K 21 , D 12 = D 21 . Furthermore, formulas (18) imply that if the packet of layers is symmetrical, then the coefficients K j j , K jl , K P j j , K P jl are equal to zero. Then, relations (17), solvable with respect to deformations, are used: The obtained values of ε i j in (19) are substituted into (17) to get To derive equations in the mixed form, the first term of (16) can be recast to its following equivalent form: ε i j from (3a) is substituted into to the first integral on the right-hand side of (16), whereas the remaining integrals take into account expressions from (19). Furthermore, M i j is taken from (20) and H i j and ε i3 -from (3c). As a result, the variational principle takes the following form Carrying out integration by parts in (21), the following variational equation, which yields a system of differential equations and boundary conditions, is obtained: Above, the following well-known differential operators were introduced: Formula (22) allows one to include relevant boundary condition, such as The system of equations governing the MM2 mathematical model of the flexible multi-layer shells was constructed for the packet of layers with the help of the Timoshenko-type hypotheses, and it looks as follows In the case of a symmetric packet, the coefficients B i j , b i j , B 31 , b 31 are equal to zero.

Modification of the Timoshenko MM2 and the ε-regularisation MM3
The Timoshenko mathematical model MM1 consists of differential equations of the 10th order and 5 boundary conditions. In other mathematical models (for instance, the Kirchhoff-Love model of the first approximation, the Pelekh-Sheremetev-Reddy-Levinson model of the third approximation), a biharmonic operator is present in the equation considering deflection. The operator improves stability and convergence of the numerical algorithm in geometrically nonlinear problems as compared to the mathematical Timoshenko model of the second-order approximation. However, construction of the computational schemes for these models does not belong to easy tasks and is more difficult than in the case of MM1. This observation is the motivation to employ the ε-regularisation [48] to establish one more variant of the improved mathematical model, called further MM3, which can be treated as a bridge between the Kirchhoff-Love (MM1) and Timoshenko (MM2) models.
In one of its variants, MM3 is a system of differential equations, in which the first equation includes biharmonic operator, -Laplace operator, ρ-radius of curvature of the contour ∂ , ν =const > 0 [49]. In addition to the problem (23), (24), one can define the problem of the form with the boundary conditions At first, the "helping problem" was considered, where the biharmonic operator was replaced with an arbitrary positive definite operator T with the natural-type boundary conditions, the energetic space of which is The fact that the operator T can be chosen arbitrarily allows for "construction" of the properties of a total algebraic system aimed at increasing the computational efficiency of the employed algorithm.
Introduction of the "helping problem" yielded a strong convergence of some series of the approximate solutions W n ε to the exact solution W 0 , which is important while constructing the practical realisation of the computational algorithm. The similar approach can be employed to the searched functions γ 1 , γ 2 , F if one substitutes the terms L 2 (γ 1 ), L 3 (γ 2 ), L 4 (F) in the system (23, 24) by the following ones: where for arbitrary i = 1, 2, 3, ε i > 0, T i serve as positive operators with natural boundary conditions, the energetic space of which is compact and Obviously, choosing the operators T, T i properly, one can achieve any required convergence of the sequence of approximations of the solutions tending to the exact value without adding any additional constraints on the input data of the considered problems.
The numerical implementation of the MM can be illustrated and discussed on the example of nonlinear differential equations governing dynamical behaviour of a flexible shallow multi-layer shell, each layer of which is made from a non-homogeneous orthotropic material within the Timoshenko kinematic model. To get a reliable approximate solution, the Bubnov-Galerkin method (BGM) and the finite difference method (FDM) can be employed.
The input problem is formulated in the following way. The aim is to find a solution to the system (27) in the space understood as a bounded part of the with the boundary ∂ satisfying the condition of the Sobolev embedding theorem [11]: and satisfying the following boundary conditions where: are known functions of stiffness [5], satisfying the following conditions A vector u = (u 3 , γ 1 , γ 2 , F) ∈ H 1 is called the generalised solution to the mentioned problem, and the following integral relations are satisfied: Besides the considered problem (25), (26), the following "helping problem" is introduced: including the boundary conditions where: 2 -biharmonic operator, -Laplace operator, ρ-radius of the contour ∂ curvature, ν=const > 0 [43]. By the general solution, one should understand a vector u ε = u 3ε , γ xε , γ yε , F ε ∈ H 2 satisfying the following integral identities: Following [44], it is not difficult to show that the functional (33) induces a nonlinear weakly-continuous operator B: H 2 → H 2 , and the functional (34) defines the linear monotonous half-continuous operator A: then for an arbitrary ε > 0 : (1) there is at least one vector u 0 ε = u 0 3ε , γ 0 xε , γ 0 yε , F 0 ε satisfying the identity (33); (2) an approximate solution to the problems (31), (32) can be found using the BGM in the following form where The theorem conditions imply that the operator A exhibits the property of coercivity. Then, owing to the Bauer's theorem [48], the system can be solved to define the coefficients of the series (36).
The coercivity of the operator A and orthogonality of the operator B allow one to write the following a priori estimate for a set of the approximate solutions (36): The above inequalities and continuity of the functional (34) as well as weak discontinuity of the functional (33) (see [50]) imply Theorem 1. One has: Indeed, the following estimations hold √ εu 3ε The inequalities (38) are yielded by (37) through the limiting transition regarding n, which allows one to choose a certain subsequence {u ε } = u 3ε , γ xε , γ yε , F ε . The subsequence is weakly convergent in H 2 and, together with the Sobolev embedding theorem, guarantees the possibility of a limit transition with respect to ε → 0 in the integral identity (37) and yields the integral identity (33) after closing the set {χ 1 } ∈ H 2 within the norm of the space H 1 .
The space can be divided into non-intersecting parts 1 , . . . , p , such that = ∪ P i=1 i . By h the maximum diameter of the partitioned spaces is denoted.
is found.
The following two functional are considered: The difference scheme is constructed by imposing a mesh on the general solution (33) in a way similar to the already discussed. To keep the properties of orthogonality of the operator B u h ε , the nonlinear operators L (u 3 , F) and L (u 3 , u 3 ) are approximated as follows L (·, ·) = (·)ȳ y (·) Wx x + (·)x x (·)ȳ y − 1 2 (·) xy (·) xy +(·)x y (·)x y +(·) xȳ (·) xȳ + Fxȳu 3xȳ ; Stability of the obtained system of equations is implied by the coercivity of the operator A u h ε and orthogonality of the operator B u h ε , i.e. the following estimations hold under applicability conditions of Theorem 1: Owing to (42), the set of points u h ε is weakly compact in H 2 , i.e. its arbitrary subset allows for extraction of a subset weakly convergent in H 2 for h → 0. Theorem 2 Any of the mentioned limits stands for a general solution to the problems (31), (32). In order to prove this, one can transit to a limit for h → 0 in the identity (39).
It is clear that conditions of the Theorem 1 can be always satisfied when the plate K x = K y = 0 is considered. Furthermore, the analogous results hold also for other boundary conditions, for instance, when the type of the boundary conditions changes along a contour [51] (for instance, this is guaranteed for an appropriate transition of boundary conditions for the Timoshenko-type model).
Remark Modifications of the Timoshenko mathematical model of the second approximation (MM2) are possible. The MM2 is created by means of employment of the kinematic hypothesis regarding tangential displacements, which yields f (x 3 ) ≡ 1 as a result of using the shell theory. However, one can impose an additional static hypothesis onto the tangential stresses σ 13 , σ 23 and take into account their variations along the thickness of each packet of the layers in the form of a squared parabola, i.e. when f ( . Such a mathematical model will be further called MM4.

Numerical investigation of stability of a multi-layer orthotropic shallow shell in the frame of various improved theories
The static stability problem is investigated similarly as described in [52], where the problem of stability was considered in terms of "load-deflection" coordinates measured at the shell centre q [u 3 (0.5; 0.5)]. The so far introduced mathematical models have been solved using the variational-difference method. The advantages of this over other possible approaches are briefly addressed below.

Minimal preliminary work (for instance, in contrast
to FEM) is required to prepare models for computer simulations. 2. Universality in terms of using various boundary conditions. 3. Possibility of development of a unique/one program devoted to investigate stability of multi-layer shells taking into account different MM. 4. Coincidence of approximations of derivatives of an arbitrary order in the given points of the space boundaries with the approximation in the internal points under imposing a mesh on the integral form of the equations of the shell equilibrium.
The following notation for the Sobolev spaces [53] is employed: (·, ·) M -scalar product, |·| M -norm in the Hibert space M. The norms in H 1 , H 2 , H 3 are introduced in the following way: 1 • . The variational-difference scheme for the Timoshenko second-order approximation model (MM2). The vector u = (u 3 , γ 1 , γ 2 , F) ∈ H 1 serves as a generalised solution to the problems (23), (24) where the identity: is satisfied for an arbitrary vector v = (ϕ, ψ, ζ, ξ ) ∈ H 1 . The functional R (u, v) can be recast to its counterpart form: where The integral identity (43) imposed on the rectangular meshω is approximated by the following total identity [54]: The difference analogue of the integration by parts [53,54] is employed: and the difference relations regarding differentiation of the mesh functions are as follows Then, the shell equations and natural boundary conditions formulated in the difference form can be written as Here, the terms of approximation of nonlinear operators corresponding to the limiting values of the variable x 1 and x 2 were not shown owing to their complexity and due to fact that the aim of this work is not to study all boundary conditions. Furthermore, the construction of MM with the boundary conditions (23a) and (23b) does not need nonlinear operators to be defined on the boundary.
The following important steps of constructing a computational scheme should be emphasised. It should be mentioned that the total identity (46) does not allow to get mesh approximation of the derivatives ∂ 3 γ 1 ∂ x 3 1 and ∂ 3 γ 2 ∂ x 3 2 while approaching the boundary layer for the boundary condition of clamping in the case of a nonsymmetric packet of layers.
In order to solve the problem, the additional conditions regarding smoothness of the functions γ 1 (x 1 , x 2 ) and γ 2 (x 1 , x 2 ) should be imposed. Namely, the values of (γ 1 ) −1, j and (γ 2 ) i,−1 in the contour-out points are yielded by difference approximation of two differential equations of the equilibrium with respect to γ 1 and γ 2 after their additional differentiation.
The differential equations are given in nondimensional forms. The relation between dimensional and non-dimensional (with bars) quantities is as follows: Note that bars over non-dimensional quantities were already omitted in equations (24). Let us consider thin shallow shells of the thickness from the interval: where R denotes the radius of the shell curvature, and the Reissner's criterion (< ρ = 60 0 ) is satisfied [8].
Recall that λ 1 , λ 2 stand for the geometric parameter characterising a ratio of the shell length and width to the shell thickness. Obviously, for each fixed value of the shell curvature, the change in λ 1 , λ 2 from minimum to maximum implies the shell thickness decrease, and hence there is no need to take angles of rotation and twisting of the normal into account. In this case, deformations of transverse shear follow: It means that transverse tangential stresses σ 13 ≈ σ 23 ≈ 0, and consequently, also transversal forces Q 1 , Q 2 ≈ 0, i.e. beginning with the defined λ 1 = λ 2 , one transits into the space of application of the classical Kirchhoff-Love model. Employing the variational-difference method and introducing the parametersλ 1 =Ā 1313 · λ 2 1 andλ 2 = A 2323 · λ 2 2 , he matrix of the system of algebraic equations is defined, where also both shell geometry and materials of layers are included by means of the ratio G i3 /ϕ 10 (i = 1, 2). Now, if λ 1 , λ 2 are increased to large values (within the possible interval), then the matrix of algebraic equations corresponding to MM1 becomes badly conditioned and convergence of the iterational process becomes too slow to employ MM1 for such kind of problems. It is clear that for large values of λ 1 , λ 2 ,one can decreaseλ 1 ,λ 2 by choosing properly the material of layers i.e. every time MM1 is employed, there is a need to follow the values ofλ 1 ,λ 2 . This drawback is cancelled in the case of using MM2, which works well in the case of both large and small values ofλ 1 ,λ 2 . This observation exhibits an advantage and universality of MM2 in comparison with the Timoshenko-type model To avoid problems with finding a solution to the system of nonlinear algebraic equations of large and badly conditioned matrix (for large values of λ 1 , λ 2 ), the iteration was employed instead of the direct Gauss-Seidel method, the main advantage of which is the selfimprovement property.
The materials of layers considered in the present study are presented in Table 1 (in non-dimensional form). To get their dimensional counterparts, the values need to be multiplied by 10 −6 (kG/cm 2 ).
In Table 2, the results of investigation of mesh convergence are reported for various numbers of nodes N (1/4 of the planform) and accuracy between iterations ε for the three-layer shells of the same thickness with curvaturesK 1 =K 2 = 24, λ 1 = λ 2 = 55.56. The exact value of the load isq = 3.1416 for the boundary condition (23a). Convergence of the numerical algorithm aimed at analysing the packet of layers of non-symmetric structure, where mesh approximations were used for differential equations with respect to γ 1 , γ 2 in the approximation of the out-of-contour nodes for the boundary condition (23a), was verified on the model problems for shells of the same structure as in the case of the symmetric packet. However, the coordinate surface was moved from the packet centre to an arbitrary point x 3 located along the shell thickness ( = h/2). Values of loads for both symmetric and non-symmetric packets of layers coincide, which validates the use of the numerical algorithm employed to the non-symmetric multi-layer shell.
3. Reliability of the numerical algorithm for MM2 with f (x 3 ) = 1 and for MM4 with f (x 3 ) = 1 − x 3 h 2 was checked by estimating deflection of the shell cen- tre u 3 (0; 0) = W 00 and comparing it with the exact Brucker's solution [55] for the counterpart geometrically linear variant of the PDEs.
One can conclude that in spite of computations of the "load-deflection" curve using different improved models, the qualitative character of the numerically obtained curves is similar. Based on extensive and numerous numerical experiments, convergence of all numerical algorithms investigating the improved MMs was found. The number of partition points N of spatial coordinates was changed with respect to the values of the curvature parameters, the geometric parameters λ 1 , λ 2 , the number of layers and their localisation regarding the reference surface x 3 = 0 in order to keep the same accuracy while comparing the numerically estimated values of critical loads.
The carried-out qualitative investigations aimed at finding the critical loads in relation to the initial data were carried out for the following numbers of mesh nodes: N × N = 11 × 11, 13 × 13, 15 × 15. One can conclude that real shells should be computed choosing an optimal mesh for particular input data, keeping the same required interval of the computational accuracy.

Comparison of stability curves of the multi-layer symmetric shells in the "load-deflection" coordinates
Since the shells composed of interlacing layers of aluminium and an orthotropic material were studied, this section presents the investigation of one-layer shells made only of either aluminium or an orthotropic material. Figures 3, 4 report the stability curves for shells made of aluminium and the orthotropic material, respectively, for fixed parameters K 1 = K 2 = 9, the orthotropic material, the results yielded by the εregularisation model (MM3) strongly differ from the results obtained by MM2 and MM4 models in "the linear part" of the q(u 3 ) curves. However, for large deflections, the difference between the results starts to decrease. In the case of one-layer shells made of aluminium (Fig. 4), the largest differences between results for MM2, MM4 and MM3 models are exhibited in the nonlinear part of the q(u 3 ) dependencies. Furthermore, Figs. 5, 6 show the functions q(u 3 ) for the three-layer shell with layers of equal thickness h = 0.046. Figure 5 presents the layers sequence aluminium-orthotropic material-aluminium, whereas the sequence orthotropic material-aluminiumorthotropic material is presented in Fig. 6. The Timoshenko model with a parabola (MM4) is very sensitive to the order of layers. In the first case (Fig. 5), the MM4 model yields the results close to the results of the ε-regularisation model (MM3). In the second case (Fig. 6), the obtained results are similar to those given by the MM2 model. It should be emphasised that the results obtained using the MM2 model qualitatively differ from the results obtained by other models. The comparison of the obtained numerical results for thinner shells of the thickness h = 0.018 and for different shell models is presented in Fig. 7 (the order of layers: aluminium-orthotropic material-aluminium).
Comparing the results shown in Figs. 5 and 7, one can conclude that the change of the layers thickness (and thus the change of parameters λ 1 , λ 2 ) yields divergence of the results given by the Timoshenko model  Numerical investigation of stability of multi-layer shells was also carried out based on simultaneous observations of stress-strain states, i.e. the stresses σ 1 1, σ 1 2, σ 1 3 for the shell deflection in the neighbourhood of the critical loads. Figure 8 present the results of the following problems: three-layer shell composed of aluminium-orthotropic material-aluminium for the fixed curvatures K 1 = K 2 = 9 and the same layer thickness h = 0.046. The results shown in Fig. 8a allow one to conclude that normal stresses σ 1 1 are similar to each other. The stresses σ 1 2 (Fig. 8b) exhibit visible differences on the shell edge and the shell centre. The stresses σ 1 3 (Fig. 8c)  All the carried-out investigations validated the statement that fundamental stresses of all tested models are uniquely defined with an accuracy allowed by the employed model for all physical-geometrical parameters with and without temperature fields. However, the shear stresses are non-uniquely defined for all investigated models. It means that the satisfaction of static conditions of the layers should be supplemented by using the discrete models while studying the stressstrain states of multi-layer shells.

Concluding remarks
The carried-out investigations yielded the following conclusions: 1. The mathematical models of geometrically nonlinear multi-layer rectangular shallow shells were derived. The governing equations were obtained in the mixed form, i.e. with respect to the stress function and deflection function. All presented models are based on the Timoshenko hypothesis. In addition, a new mathematical model with εregularisation was constructed, and the theorem on existence of its general solution was formulated and proved. 2. The conservative difference schemes were developed for the boundary value problems, based on the fundamental variation equations of the shallow shells. The method based on the difference approximation of differential operators of the odd order yielded by the theory of multi-layer shells with nonsymmetric packet of layers was developed. 3. Investigations devoted to feasibility of the constructed mathematical models were carried out. Reliability of the obtained results was analysed and validated by means of comparing the obtained results with the analytical and/or numerical results obtained by other authors. 4. The carried-out numerical experiments showed convergence of numerical algorithms of all developed mathematical models. It was demonstrated that the stability curves in the "load-deflection" coordinates depend on the layers position, layers thickness and the layers material. 5. The one-layer shells made from two different materials were investigated within all constructed models. It was detected that the mathematical models of the Timoshenko type with f ( shows that normal stresses yielded by all models are very similar to each other. Also, the stresses σ 1 2 are qualitatively similar for all models, but the ε-regularisation model exhibits a difference on both the shell edges and the shell centre. The largest difference was obtained for the stress σ 1 3, where for MM2 (MM4) and MM3, the detected difference is both qualitative and quantitative.