Contact interaction of two rectangular plates made from different materials with an account of physical nonlinearity

A mathematical model of a contact interaction between two plates made from materials with different elasticity modulus is derived taking into account physical and design nonlinearities. In order to study the stress–strain state of this complex mechanical structure, the method of variational iteration has been employed allowing for reduction of partial differential equations to ordinary differential equations (ODEs). The theorem regarding convergence of this method is formulated for the class of similar-like problems. The convergence of the proposed iterational procedure used for obtaining a solution to contact problems of two plates is proved. In the studied case, the physical nonlinearity is introduced with the help of variable parameters associated with plate stiffness. The work is supplemented with a few numerical examples. Both Fourier and Morlet power spectra are employed to detect and analyse regular and chaotic vibrations of two interacting plates.


Introduction
Many researchers have indicated that several structural materials may exhibit different tensile and compression response under flexural testing. For instance, composites fabricated from carbon fibres pre-impregnated by epoxy polymer usually present different elastic responses under tension and compression in the principal material directions. This phenomenon has been observed for the first time by Jones [1], who introduced the name of multi-modulus materials.
The multi-modulus materials, such as concrete, are widely employed in civil engineering [2]. While investigating radial and tangential stresses around boreholes, it was found that a slope of the linearized relationship in compression is almost always higher than that in tension [3]. Different stress-strain characteristics of cast iron and Poisson's coefficient in tension and compression were reported by Gilbert [4].
Patel et al. [5] carried out both theoretical and experimental investigations of material properties of a calendared ply of numerous composites widely used in the body and belt of radial tires. They detected bi-modulus behaviour under tension and compression.
Barak et al. [6] employed electronic speckle pattern interferometry to measure concurrently the compressive and tensile strains in cortical bone beams tested in bending. It was reported that the tensile Young's modulus is slightly (6%) greater than the compressive Young's modulus.
Bertoldi et al. [7] developed a micromechanical model for nacre (mother of pearl) by using a homogenization approach. In the small-strain regime, the nacre exhibited strong difference in directional stiffness and bi-modular effects (different Young's moduli in compression and tension).
It should be emphasized that even graphene, the strongest material having minimum electric resistance, shows greater compression modulus than tension modulus [8,9]. Graphite, which is a descendant of graphene, was investigated by Yao et al. [10] and different tensile and compressive elastic moduli were measured in this study. In addition, buckling tests were performed while studying the buckling behaviour of a graphite rod. The analytical and finite element numerical models demonstrated to be reliable.
Yao and Ye [11] proposed the analytical/numerical model of thermal stresses for a structure with different moduli.
In the last century, a few mathematical models were developed to study behaviour of the multi-modulus materials. The first one was proposed by Bert [12], who studied a certain composite material exhibiting different Poisson compliance when loaded transversally to the fibres than when loaded parallelly to the fibres. The proposed method has been commonly used in multilayer composites [13][14][15][16]. The second widely accepted model was proposed by Amburtsmian [17] and was validated for isotropic materials. It allows for estimation of different tension and compression moduli based on the positive/negative sign of the fundamental stresses, which play an important role in engineering application.
Furthermore, Jones [18] proposed the modified and improved model called the weighted compliance matrix (WCM) material model, which satisfies the criteria for isotropic and orthotropic bodies under plane stresses. The state of the art of the multi-modulus materials was presented, for instance, in reference [19] and thus is omitted here.
The contact problems have a long history in mechanics, and at the beginning, mathematical approaches were applied to solve them. The problem of a point contact between elastic bodies was formulated for the first time by Hertz [20]. Nowadays, a rapid development of technology and fabrication tools has attracted numerous researchers to revisit modelling of the contact problem, which stands for one of the most challenging problems in the field of mechanics of deformable rigid bodies. The complexity of the problem requires rather sophisticated mathematical approaches aimed at finding precise, reliable, and validated solutions. In particular, the complexity and difficulty of the mentioned problems increase while dealing with the contact interactions of structural members made from materials having different moduli.
Duvant and Lion [21] proposed the first attempt to study contact problems with the help of evolutionary variational inequalities. A deeper state of the art of mathematical approaches can be found in references [22][23][24]. More recently, quasi-static and dynamical frictional non-standard contact states belonging to a new class, including also effects of memory, have been studied in piezoelectric materials [25,26].
The so far briefly reviewed state of the art emphasizes the importance of analysing the phenomena associated with contact problems of plates taking into account the different moduli of the used materials. Our work brings a contribution to this problem by derivation of a mathematical model of the contact interaction of two plates, where the physical nonlinearity is also taken into account. The paper is organized in the following manner. Section 2 of the paper deals with derivation of the partial differential equations (PDEs) governing the dynamics of two plates made from bi-modulus material. In addition, the plates are of variable thicknesses and are separated by a clearance. Reduction of PDEs to ODEs is described in Sect. 3, including a numerical example. Section 4 presents theorems associated with the MVI (method of variational iterations) convergence putting emphasis on our formulated theorems. Contact interaction of two square plates and numerical examples is reported in Sect. 5. Section 6 presents dynamic problems of contact interaction with employment of Fourier and Morlet wavelet spectra aimed at investigating regular and chaotic vibration of the structural package. The paper ends with a Sect. 7 containing concluding remarks associated with the obtained results.

Statement of the problem
We consider a problem dealing with one-constrained mechanical interaction between two rectangular plates of variable thickness. The plates are treated as thin, and their stress-strain state (SSS) can be described by the classical Kirchhoff theory supplemented with physical nonlinearity and considered in the frame of the theory of small elastic-plastic deformations. We assume that the contact pressure (normal to the surface of generated stress) is much less in comparison with the normal stress occurred in the plates cross sections, and hence, a free slip of the contacting plates is allowed. A choice of the classical theory of plates (the Kirchhoff model) is motivated by an observation that, in comparison with the transversal crushing in a contact zone, the deformation of the transverse shear has a negligible influence on both the SSS and the distribution of the contact pressure [27]. This validated observation stands for us as the main purpose of developing and further employing our theoretical background.
Let us present a system of the input PDEs for two contacting plates in the following form [27]: where w i (x, y)-vector functions; q i (x, y)-vector of the external load; i = 1, 2-number of a plate coinciding with a positive sense of the normal; , (x, y) ∈ Ω i are the functions governing forms of external surfaces of the plates, which allow to take into account changes in thickness of each plate (Fig. 1). The influence of different plate moduli on the plate behaviour is investigated based on Zhao and Ye [28] as well as He et al. [29] proposals. Each plate is divided into two subspaces by means of employing a concept of the so-called elastic neutral layer z = c i (x, y):  [29], is defined by non-constant elasticity modulus and Poisson's coefficient, and the following relations hold The following notation is employed: z = a i (x, y) and/or z = b i (x, y), (x, y) Ω i denote equations of external surfaces of the plates, which allow for inclusion of variable thickness of each plate; i ) describe different stiffness modulus and Poisson's coefficients, respectively; e (i) i stands for intensity of deformation of i-th plate.
It should be emphasized that the method of variable parameters of stiffness [30], employed in our physically nonlinear problem, is based on observation i ) are non-constant parameters dependent on the deformation state in a plate point, and they are governed by the following formulas where K 1i = const. The theory of deformation yields the following shear modulus where σ (i) i , e (i) i stand for intensity of stresses and deformations of a plate, respectively, and The value of e (i) zz is defined through the condition of a plane stress state (s zz = 0): Two types of boundary conditions are analysed, i.e.
(i) simple support (ii) rigid clamping Note that in mathematical models (1), (4), (9), (10) and (1), (5), (9), (10), the contact pressure is proportional to the transversal crushing w 1 − h 1 − w 2 in a zone of contact between plates, i.e. we have and the function ψ has the following form Here h 1 stands for a clearance between plates, and k is the stiffness coefficient of the transverse crushing of the plates. Due to occurrence of the multiplier ψ in (1), the system of equations becomes nonlinear, and consequently the entire problem becomes nonlinear. The design nonlinearity governs the possibility of switching on/off of the one-sided constraint. Furthermore, inclusion of the physical nonlinearity in the considered problem makes it impossible to solve the problem via classical analytical, semi-analytical, or numerical approaches. Therefore, we employ here the method of successive approximations matched with a simultaneous process of improvements of the borders of the contact zones.
Observe that formula (11) holds for the case of a contact between plates having the same thickness h and the same parameters E and k. In the case of contact problems of Kirchhoff's plates, Winkler-type coupling between the crushing zone and the contact pressure is utilized. If there is a thin layer between plates, it can be also taken into account by means of changing the parameter k.
If the initial plates configuration (function h 1 ) and the employed load correspond to lack of contact between the deformed plates, then ψ ≡ 0, and the system (1) is split into two separate subsystems. Otherwise, two subsystems (plates) are coupled. Substituting (11) into (1) yields Observe that the system (13) is of the eight orders and hence it should be studied together with the boundary conditions (9), (10). It includes both design and physical nonlinearity. In order to solve the derived design-nonlinear problem (13), one has to construct an iterational process allowing (for each step of loading) for solving only one of the equations of the system (13). This procedure yields a decrease in the systems order twice for a twolayer package and, consequently, into n times while dealing with a package of n layers. The iterational procedure has the following form The system (14) should be supplemented with the corresponding boundary conditions (6) or (7) for the i-th plate. If the values of E and h are different, then for the upper/lower plate we use notation E 1 , h 1 /E 2 , h 2 , and we have In the case of fixed/frozen contact zone, Eq. (14) can be solved using one of the methods devoted to reduction of PDEs to ODEs. Those methods are briefly discussed in the next sections of the paper.

Reduction of PDEs to ODEs
Nowadays, there exist a few methods aiming at reducing the PDEs to ODEs. The principal concepts and ideas of those methods will be discussed using the operatortype equation of the form where A-operator of a boundary problem including both differential operator and boundary conditions; q(x, y)-known function (vector function of the system of equations); w(x, y)-function being searched (vector function of the system of equations).

The method of Kantorovich-Vlasov (MKV)
Owing to the MKV, the following form of solution to (16) can be assumed [31,32]: where ϕ i (x) are given functions satisfying the boundary conditions (9) or (10), whereas ψ i (y) are the being searched functions defined by the following system of projected equations Here (Au N − q, ϕ j (x)) stands for notation of a scalar product in the space is the space in which the exact solution w 0 to Eq. (16) exists. The procedure (18) reduces the problem of PDEs to ODEs by employing the Bubnov-Galerkin method regarding one of the coordinates. This procedure, in fact, generalizes the classical Bubnov-Galerkin method [33].
Disadvantages If the expected solution (17) should obey the same symmetry property with regard to all variables, then the MKV method sometimes cannot exhibit this property due to the limited number of terms used in the series. In other words, in case of taking a limited number of the series terms, the mentioned symmetry can be lost.
Advantages Independence of the basic functions versus one of the variables (here y) on the form of the boundary conditions.

Method of Vaindiner (MV)
In this case a solution to Eq. (15) is assumed to be as follows [34]: where ψ i (y) and u i (x) are given functions, ϕ i (x) and v i (x) are the being searched functions satisfying the system of projected equations of the following form where Drawbacks There is a need to choose a full coordinate system of equations; there exists a dependence of accuracy of the obtained solutions on a number of equations in the system of approximating equations 2N , i.e. the number of equations is twice increased in comparison with the MKV.
Advantages In spite of the same advantages as those exhibited by MKV, there is a possibility of obtaining a symmetric solution with respect to all variables, if the being searched exact solutions have the mentioned symmetry.

Method of variational iteration (MVI)
This method presents a modification of the MKV [35]. Owing to MVI, a solution to Eq. (16) can be searched in the form of (17), where the functions ϕ i (x) and ψ i (y) are defined through the following equations (22) in the following way. A certain system of N functions with regard to one of the variables, let us say ϕ 0 j (x) j = 1, 2, . . . , N , is given. Then, the system (21) defines a system of functions ψ (1) j (y), which are substituted into the system (22), and a new choice with regard to the variable x is defined as ϕ (2) j (x). Then, using the chosen functions, the system (21) yields a new function ψ (3) j (y) with respect toy, and so on. Truncating the process of searching for the functions ϕ j (x) and ψ j (y) on the k-th step, we define the following function which serves as an approximate solution to Eq. (14). Disadvantages The accuracy of the obtained approximate solution depends on a number of the series terms of (17).
Advantages There is no need to define the initial approximations satisfying, for instance, the boundary conditions; symmetry of the approximate solutions, if it is a priori known, holds; it is possible to achieve the maximum allowed accuracy of the approximate solutions even for a bounded number of terms approximating the solutions.

Method of Arganovskiy-Baglay-Smirnov (MABS)
The so far discussed method of variational iterations (MVI) coincides with the first term of the series of the method developed by Arganovskiy and Baglay [36], Baglay and Smirnov [37]. The formal scheme of the MABS can be briefly summarized in the following way. A solution to Eq. (16) in the first approximation (N = 1) is searched in a way similar to the MVI, i.e. we take The new equation is defined as follows i.e. we have changed the right-hand side of Eq. (16). Equation (16) is solved again with the help of MVI, and its first approximation yields The next new equation follows: Aw 3 (x, y) = q(x, y)− Aw 1 (x, y)− Aw 2 (x, y), and then one employs the MVI again in the first approximation, and so on. Finally, the following series is used as the input solution: Drawbacks This method cannot be used to solve a class of problems (for instance, if the boundary conditions are changed in a rectangular plate).
Advantages The advantages are analogous to those of MVI. Furthermore, each of the defined right-hand sides of the type of Eq. (25) can be solved by the MVI only in the first approximation, which generates algebraic systems of a small dimension.

Combined method (MC)
This method is a combination of the MV, MABS, and MVI. Owing to the MC, the approximate solution is searched in the following form where ψ 0 (y) and u 0 (x) are assumed as known functions; ϕ (1) (x) and v (1) (y) are the being searched functions defined by the projected Eq. (16). However, on the contrary to the MV, where the searched functions ϕ (1) (x) and v (1) (y) have been treated as the final ones, here they are taken as known functions used to define the new ψ (2) (y) and u (2) (x) unless the accuracy between u (k) (x, y) and u (k−1) (x, y) exceeds the assumed magnitude. We are looking for a solution of the equations defined through the so far described way, i.e.
If the error between u 1 and (u 1 +u 2 ) is less than a given quantity, then the process is terminated. If not, then the following equation is constructed and w 3 (x, y) is defined, and so on. The final solution to Eq. (16) is defined by the following series Advantages This method includes all advantages of MV and MABS; furthermore, it may be used for a wide class of problems (for example, for the problems of rectangular variable boundary conditions along their contours). The key idea of this matching is the application of the procedure of MABS to the MKV without the employment of the method of variational iterations (MVI).
Owing to the method, an approximate solution of Eq. (16) is found based on the Kantorovich-Vlasov method, assuming its following form where ϕ i (x) are the given functions satisfying the corresponding boundary conditions (9) or (10), ψ i (y) are the searched functions defined by the system of projection Eq. (18). However, if the functions w 1 (x, y) found in the MKV are the final ones, they can be employed to construct the operator Aw 1 (x, y) on the right-hand side of equation which is solved using the MKV. If the error between w 1 and w 1 + w 2 is smaller than a given quantity, then the iterative process is stopped. If this is not the case, than we construct the following equation (34) and we define w 3 (x, y) with the help of the MKV, and so on. The terminal/final solution to Eq. (16) is given by the following formula 3.7 Matching of the methods of Vaindiner and the Arganovskiy-Baglay-Smirnov (MV + MABS) In this case the method of Vaindiner (MV) is supplemented with the procedure of the MABS but without a procedure of MVI. A solution to Eq. (15) is searched in the following form: where ψ i (y) and u i (x) are the given functions, ϕ i (x) and v i (x) are the searched functions defined by the system of projection Eq. (18). However, while the found function w 1 (x, y) is treated as the final one in the MV, it is used to define the operator Aw 1 (x, y) on the right-hand side of the following equation which is solved by the MV. If the error between w 1 and w 1 + w 2 is smaller than a given quantity, then the iterative process is stopped. Otherwise, we construct the following equations (38) and we define w 3 (x, y) with the help of the MKV. As the final solution of Eq. (16), we take the following 3.8 Matching of the methods of Vaindiner with the method of variational iterations (MV + MVI) In this case the method of Vaindiner (MV) is extended to include a procedure of the variational iterations. The approximate solution to the problem is searched in the following form where ψ 0 (y) and u 0 (x) are known functions, ϕ (1) (x) and v (1) (y) are the searched functions defined by the system of projection Eq. (16). However, if the functions ϕ (1) (x) and v (1) (y) found with the help of MV are the final ones, they are taken as known functions, and the being searched solution takes the form w (2) (1) (y). (41)  Employment of the MVI yields new functions ψ (2) (y) and u (2) (x), and so on. The iterative procedure is carried out until the error between w (k) (x, y) and w (k−1) (x, y) is smaller than a priori given value.

Numerical example
In order to compare the accuracy of each of the discussed methods, we consider the problem of a linear bending of the Germain-Lagrange plate governed by the following PDE which has known exact solution.
The differential operator (42) and the employed boundary conditions are presented in their counterpart nondimensional form x =xa, y =ȳb, w =wh, and the bars over the non-dimensional quantities are further omitted. Note that a system of integral-differential equations yielded by the Kantorovich-Vlasov procedure can be solved through the method of finite differences (FDM) with the successive employment of the Gauss method while solving a system of algebraic equations produced during looking for the solutions. The solution to the stated problems has been solved using unitary and doubled accuracy. The interval [0, 1] has been divided into 30 nodes. Further increase in nodes number has not changed the obtained results. The relative error of computations has been traced using the formula The obtained results are reported in Table 1 forq(x, y) = 50.
The so far described methods as well as the results of the numerical experiment allow one to conclude that the MC has the highest accuracy in comparison with other methods and it yields the results being close to the exact solution. However, a drawback of this method is associated with the increase in the number of differential equations, and hence, one more iteration should be involved to get reliable results, which implies an increase in the computational time. However, the method of variational iterations is free from this drawback, and the yielded results are also close to the exact ones. Therefore, particularly in the case of complex problems, the MVI is recommended.

Theorems on convergence of MVI
In the case of a fixed contact zone, the procedure of the MVI [35] is constructed (see the comments given in Sect. 3.3) together with the method of variable parameters of stiffness (MVPS) [30] unless the required accuracy is achieved. A successive improvement of a contact zone is achieved due to the method of direct iteration. Then the solution procedure is repeated, i.e. we employ three iterational procedures nested in each other. Scheme of MVI can be formally described in the following way. We are aimed at finding a solution to equation where A stands for a certain operator defined on the manifold D(A) of the Hilbert space L 2 (Ω); q(x, y) stands for a given function of two variables x, y; w(x, y) is a searched function; Ω(x, y) is a space associated with variations of x and y.
If Ω(x, y) = X × Y (X-a certain bounded set of variables x; Y -a bounded set of y), then a solution to Eq. (46) has the following form where the functions u i (x) and v i (y) are defined by the following system of equations in the following way. A certain system composed of N functions with respect to one of the variables, for instance, u 0 1 (x), u 0 2 (x), . . . , u 0 N (x) is given. Then, the first N equations of the system (48) Next, the obtained functions are employed to create a new set of functions x − u 2 1 (x), u 2 2 (x), . . . , u 2 N (x), which is further used to construct a set of new functions with respect to the variable y, i.e. v 3 1 (x), v 3 2 (x), . . . , v 3 N (x), and so on. In the case of the iterational procedure MVI [35], proofs of the theorems constituting the theoretical background of the MVI convergence were given for the problems of the theory of plates.

Theorem 2 Let each element of the basis system of the
where is a basis system in the space The method of variational iterations (MVI) removes the necessity of constructing the system of approximating functions, on the contrary to the widely employed Bubnov-Galerkin and Ritz methods [38,39]. The initially given arbitrary functions (perhaps satisfying the well-known conditions of smoothness) are modified in the process of computations through MVI beginning with a solution of the input system of PDEs governing the plates behaviour. Therefore, in a natural way, the MVI constructs the approximating functions in the form of a product of two functions, where one of them depends on one argument in a way similar to the Fourier method. This implies the necessity of solving two ODEs instead of one. In fact, the method of variational iterations generalizes the Fourier method.
The most fascinating result is that the procedure of the MVI is converged even if functions not satisfying the boundary condition (9), (10) in the form: w 1 (x) = 1; w 2 (y) = 1 are taken as the initial approximations.
The MVI procedure was stopped already on the fourth step of the variational iterations. Remarkably, almost accurate coincidence with the exact solution has been obtained, i.e. the boundary conditions (6)/(7) correspond to the difference between the approximate and exact solution in percentage of 0.1/0.6%.

Convergence theorem
In what follows, we formulate a theorem (without a proof) regarding convergence of the iterational procedure (14) associated with solving our formulated problem aimed at designing a nonlinear system composed of two plates made from materials of different moduli.
Let R 2 be an Euclidean plane with a Descartes basis: Ω i R 2 -area on a plane with the boundary We consider a space plate-type object composed of two contacting plates (Fig. 2) obeying Kirchhoff's hypotheses and governed by the following PDEs: with the boundary conditions (10). The function defining the contact zone Ω * of the plates is as follows and q 1 (x, y), q 2 (x, y) are the functions of external load acting on the first and the second plate, respectively; the operators A i (w i ) have the form (5), which in the case of elastic problems is governed by (4); · Anorm in a normalized space A; (·, ·)-scalar product in the Hilbert space L 2 . The employed notation of the functional spaces coincides with that given in reference [40], whereas E 1 , E 2 stand for different moduli in tension and compression. In order to solve the problem (51), (10), the following iterational algorithm is employed: Then, the following theorem holds. Theorem Let Ω i , (i = 1, 2) are bounded spaces and let their boundaries ∂Ω i satisfy conditions of the Sobolev embedding theorems [40], where Ω * stands for a measurable space, q i (x, y) ⊆ L 2 (Ω i ) and besides, there exist real constants c i > 0, D i > 0 such that

Remark 1
The method of the variable stiffness parameters allows one to study SSS of the plates and plate-type

Remark 2
The given theorem implies the theorem reported in reference [39] valid for a material with one modulus.

Contact interaction of two square plates
In this section we consider bending of two squared plates made from materials having different moduli in tension and compression, taking the physical nonlinearity into account (see Fig. 2). If values of the material parameters E + , E − , ν + , ν − are slightly changed, the change in the elastic neutral layer (surfaces of the zeroth order) can be neglected [28,29]. This observation has been taken into account in our further considered computational examples.
In order to investigate the stress-strain state of the plate-type structures made from materials of different moduli, one needs to use two different dependencies of the stress intensity on deformation. Let σ ) be dependencies of the stress on deformation for the i-th plate. ("+" refers to tension, whereas "-" refers to compression.) Depending on a sign of the averaged stress of i-th plate σ i ) corresponding to either tension or compression for each of the plates. In this case, the previously implemented algorithm, including the method of variable elasticity parameters, method of variational iterations, iterational method of variational iterations, and iterational method of contact interaction, is employed.

Computational examples
We consider an example of a computation of the contact interaction between the square plates made from materials of different tension/compression moduli.
The plate material is a pure aluminium, and its nonlinear diagram is governed by the following formula where e i S = σ S 3σ 0 is a conventional intensity of plasticity deformation, and hence, While approximating the σ i (e i ) dependence shown in Fig. 3, one may use the following formula The coefficient b can be identified assuming that the curve L goes through the point (σ This allows for more feasible tracing of the velocity of approximating σ s by σ * i with an increase in e i . In the case of a material with different moduli (Fig. 4) and with an account of σ i (e i ) in the case of physical nonlinearity, we have (59) The conventionally recognized notation ± for σ i , e i , e is denotes the values of the stress, deformation and plasticity deformation, respectively, for the case of tension (+) and compression (-). The coefficient α may have values larger or smaller than 1 in both cases of deformation, i.e. tension or compression.
As an example, we consider two square plates simply supported along their contours [see (10)] for λ = a b = 1 and for their relative thickness a/ h = 12.5. Each plate is made from aluminium and exhibits the stress plasticity threshold σ S = 1023 × 10 5 Pa, G 1 = G 2 = 0.383 × 10 11 Pa (here stresses are measured with respect to G 1 (h/a) 2 ), the non-dimensional deformations e = e i (a/ h) 2 , whereas the plasticity deformations are described by the formula e S = e i S (a/ h) 2 . In reference [41], it has been shown that the properties of aluminium are well approximated by the dependencies σ i (e i ) through formulas (55), (56) and they are validated experimentally. Owing to the data reported in this reference, we take k = k 1 = 2G 1 , and the σ i ) function is governed by the following formula the algorithms developed by us for estimating the dependence σ i ) can be given in an arbitrary way, i.e. either in a tabular form (digits) or in an analytical way. The volume of each plate is covered by a mesh B i with respect to x, y, z of the resolution 30 × 30 × 8. In each point of two volumes, the following quantities have been computed σ av , E i , ν i (i = 1, 2) as a result of employment of the formulas reported in Sect. 2, and hence, the iterational procedures have been constructed.
The load-deflection dependence in the centre plate of the upper plate q(w 1 (0.5, 0.5)) is reported in Fig. 5, where q(x, y) = q(x,y) The values of α − are shown in the figure, whereas α + = 1.0 for all investigated variants. Observe that an increase in the loading yields essential differences in SSS in comparison with decreasing α − . The value of the threshold load decreases on amount of 2.1% with a decrease in σ s on amount of 5%. Now we consider the SSS of the contact two-layer construction simply supported along the contour (9). Here the upper plate is made from a material of different moduli and physically nonlinear (solid curves), whereas the lower plate is made from one-modulus material (dashed curves). Relations q(w 1 (0.5, 0.5)) are reported in Fig. 6.
The value of deflection, for which a contact occurred, is denoted by a dot (w 1 = 0.01). In the case of contact problems, the differences exhibited by the curves q(w 1 ) are essentially less than for the one-layer plates for α − = 0.8; 0.95; 1.0. Analysis of the SSS of the plates implies that although the qualitative behaviour of the one-modulus (α − = 1) and different moduli Fig. 6 The load-deflection dependence of the upper plate (only upper plate is made from material of different moduli) (α − = 1) materials is the same, there is an essential quantitative difference. The ∂w ∂ x for different modulus α − = 0.8 of the first plate and one modulus of the second (lower) plate is shown in Fig. 7a, b, respectively. Figure 8a, b shows the surfacesM x (x, y) for the upper plate for the same parameters α − = 0.8, and 1.0, respectively. Analysis of the plates SSS implies that their qualitative behaviour is similar in both cases (a − = 1), and (α − = 1), but there is an essential quantitative difference.

Dynamics of a contact interaction
It should be emphasized that the phenomenon of contact interplay exhibited by mechanical systems composed of multi-layer plates coupled by boundary conditions is reported in this paper for the first time.
The further employed wavelet-based analysis allows one to study different kinds of vibration character in a given time interval. It allows also for detection of a birth/death of a frequency and time of its occurrence/vanishment. Although there exist different types of wavelets [42,43], the Morlet wavelet is the most suitable for our purpose. The analysis of advantages and disadvantages of different wavelet types to be employed in the study of shells can be found in reference [44].
A phase synchronization means phase locking of chaotic vibrational signals in time, although amplitudes of the mentioned signals remain uncoupled with each other and exhibit chaotic behaviour. Furthermore, phase locking implies frequency locking (res-onance phenomenon). The wavelet surface w(s, t 0 ) = w(s, t 0 exp[ jϕ s (t 0 )]) characterizes the system behaviour in each time interval s at an arbitrary time instant t 0 . The quantity w(s, t 0 ) describes the occurrence and intensity of the corresponding time interval s at the time instant t 0 . The integral energy distribution of the wavelet spectrum with respect to time intervals E(s) = w(s, t 0 ) 2 dt 0 is introduced. The studied phase is defined as ϕ s = arg W (s, t) for each time interval s, i.e. one can characterize behaviour of each time interval s with the help of the associated phase ϕ s (t). In particular, we have employed the Morlet wavelet while investigating the synchronization phenomenon of the considered mechanical structures.
Chaotic phase synchronization plays a key role in modern theory of nonlinear vibrations. The occurrence of the chaotic phase synchronization has been experimentally validated in radio-frequency generators, electromechanical oscillators, lasers, heart beating, and many others. It has also a significant influence on robustness of information transmission with the help of deterministic chaotic vibration.
Equations governing the dynamics of two-layer package of physically linear elastic plates have the following form (see Fig. 2) In a particular case of one-modulus (E System of Eq. (62) has been obtained by the following introduced relations between dimensional and non-dimensional parameters: where a, b-plate dimension with respect to x and y, respectively; ttime; ε-damping coefficient; w-deflection; h-plate thickness; ν = 0.3-Poisson's coefficient; g-the gravity of Earth; E-Young's modulus; q(x, y, t)transverse load; γ -specific material weight. Note that bars over non-dimensional quantities have been omitted in (62) for simplicity reasons.
As boundary conditions, simple support along a contour is taken: The initial conditions are as follows We take the control function ψ = 1 2 [1 + sign(w 1 − h * −w 2 )]; Ψ = 1, if w 1 > w 2 + h k , i.e. there is contact between plates; otherwise Ψ = 0; w 1 , w 2 -function of deflection of upper and lower plates, respectively; K -stiffness coefficient of transverse clamping in the contact zone; h * -clearance between plates, and q(t) = q 0 sin(ω 0 t).
A solution to system (60)-(62) is being searched by the Faedo-Galerkin method, and the problem is reduced to investigating nonlinear ODEs (the Cauchy problem). The function w i , i = 1, 2 is approximated by a function being a product of two functions (one depending on time and the other depending on coordinates) of the following form The derived Cauchy problem has been solved using the fourth-order Runge-Kutta method.
As an example, let us consider vibration of the plates under the sinusoidal and uniformly distributed load for the following fixed parameters: q 0 = 2, K = 5000, h * = 0.01, ε = 0.5, ω 0 = 5.0887. (The last one stands for the fundamental frequency of linear plate vibration.) It has been shown that to verify convergence of the numerical results obtained by the Faedo-Galerkin method for N = 1, 3, 5, it is sufficient to take only three first terms of the series in formula (65). It has been found that once a contact between two plates takes place, vibration become chaotic. The character of plates vibration has been investigated with the help of qualitative methods of nonlinear dynamical systems, i.e.
we have employed Fourier and wavelet-based analysis, studied contact pressure between plates, and monitored the phase difference of plate vibration.
Graphs of simultaneous vibration of the centres of both plates w(t), their contact pressure Pq(t), time his-tory of phases difference Δϕ s (t) as well as the Fourier power spectrum F(ω) and 2D/3D wavelet spectrum for each of the plate, respectively, are presented in Fig. 9.
Time histories of plate deflections (a) and time histories of contact pressures (b) exhibit complex chaotic interaction of the plates. A dark colour used in the plots showing the phase difference (c) exhibits the frequencies zone ω ∈ [4,8] where the phase synchronization has been detected. Vibrations exhibit clearly the frequency ω 0 , i.e. the vibrations take place on the fundamental frequency of free vibrations of the two-layer structural package. This is validated by the Fourier spectrum (d, g) and the Morlet spectrum (e, h, f, i). The latter exhibits energetic component of each of the frequencies at a given time instant.
Remark 1 It should be mentioned that the so far solved dynamical problems and used methodology can be successfully employed also to solve static counterpart problems. The method is known as the relaxation method and has been used by us earlier to solve different nonlinear problems [45]. From the mathematical point of view, the relaxation method belongs to a class of iterational methods devoted to solve a system of algebraic nonlinear equations where each time step creates a new approximation to a being searched solution. This method is of high accuracy, and it is robust to the choice of the initial approximation. This observation is motivated by physical correspondence of Eq. (50) to vibrations of plates embedded into a viscous environment.

Remark 2
In the case of the studied nonlinear dynamical problems, after transition from PDEs to ODEs, one can (on each time interval) improve the system of approximating functions w n (x, y), where a solution to the dynamical problem can be presented in the following form by using one of the earlier described methods (Sects. 3.2-3.5). This approach can be used particularly to the problems where the boundary conditions are changed due to vibrations.

Concluding remarks
The following general conclusions can be formulated as a result of the carried out research.
1. The proposed approach is convenient for analysis of SSS of multi-layer plates and beams with welded layers and with the occurrence of separated contacting zones due to the technological and exploitation reasons. 2. The method offers a possibility of studying the constructions with contacts between layers, when the contact condition may depend on coordinates and all non-homogeneity properties of the studied contact. 3. Important results yielded by our theoretical investigations concern removal of the function describing the contact interaction between the structural members from a number of unknowns. 4. The given algorithm allows one to study SSS of the structures, in which the dependencies σ i ) can be given either in digital forms or can be formulated analytically. 5. In the case of a two-layer plate-plate package with design nonlinearity, phase synchronization takes place in the frequency interval ω ∈ [4,8], where there is lack of synchronization (bright zones in the graphs of the phase difference) in certain time intervals, but amplitude locking phenomena appear instead.
x q 2 l 2 1 l Fig. 10 Contact between a thin layer of the thickness 2l 1 and a die the problem of contact of a layer of thickness 2l 1 with a rigid die. Let the layer be under pressure q in a zone of the length 2l symmetrically distributed with respect to axes x, z (see Fig. 10). We restrict the considerations to the plane case and we develop pressure into a trigonometric infinite series of the following form The stresses are derived in reference [46], and they take the following form We take q(x) = −σ z (x, l 1 ). In the case of a plane in a deformable state, we have It should be emphasized that the so far described model of the transverse clamping in contact zones can be also employed in frame of the classical kinematical Bernoulli-Euler hypothesis. The latter statement has been illustrated and discussed in reference [47], where both normal load and contact pressure satisfy the same requirement: they are considerably less in comparison with normal stresses in transverse shell cross sections and, consequently, they can be neglected in Hook's law relations.
Integrating (A.2) with respect to z, and taking into account that l 1 0 σ x dz = 0 and since there is no load in the x direction, we obtain Both nominator and denominator of (A.3) are developed into a series with regard to βl 1 ≡ c, and hence, (A.3) yields Now, introducing the notation 2l 1 = h, we obtain a formula for coupling of the transverse clamping Δ = 2u z (x, t) and the function q(x) of the following form x (x, t) + · · · (A.5) Next, we consider clamping of the same layer two smooth, symmetrically located with respect to axes x, z (Fig. 10), dies. In this case, a solution to the contact problem, where the function u z (x, l 1 ) is given on the interval −l ≤ x ≤ l, satisfies the following equation was reported in reference [48], where the approximation u L(u) = 1 A + ∞ i=1 B i u (2i) ( p) was used to obtain the following relation: stands for the derivative of 2i-th order (Fig. 11).
In what follows, we find the solution in an explicit form of the mentioned formulas [48]: Owing the evenness of the function K * ( p), we have K * ξ −x One can verify that the relations (A.5) and (A.12) are mutually inverse. They imply that when there are no bending deformations and elongations of the middle beam surfaces, in the two so far studied tasks, a coupling between clamping and contact pressure differs from Winkler's model only by the derivative of the fourth order multiplied by small coefficients. Even for h/l = 1 (l stands for a characteristic dimension of a contact zone), the first coefficient of the so far mentioned coefficients is equal to 1/720, where usually h << l.
The so far well-judged considerations yield a conclusion that Winkler's model (A.12) exhibits high enough accuracy. The latter conclusion has been also validated by the authors of reference [49] under support of the asymptotic analysis as well as by the author of reference [47] where a wide class of contact problems of symmetric shells has been solved.
If one takes into account tangential deformations yielded by bending and elongation of the contacting beams, the coefficient Δ changes, and small additive terms appear. For instance, the formula of clamping in the case of plates yields q k = 3(1−υ) 4(1+υ)(1−2υ) E h Δ + o(h 2 Δ 0 w) . In the case of beams, taking into account the Poisson's coefficient υ = 0, it yields q k = Analogous formulas have been also obtained and studied in references [50,51]. They differ only in the coefficient K standing by E/ h, which in all related formulas is close to one (1).
Taking into account the so far carried out considerations and owing to the mentioned dedicated references while solving the contact problems of Kirchhoff theory of shells, Winkler's model of coupling between the clamping and contact pressure has been validated.
In what follows, we express the contact pressure between the Euler-Bernoulli beam and the dies by the difference Δ = w − a, where w stands for a normal displacement of the middle surface [49,[52][53][54], and the following simple formula holds ) , a > 0. (A. 13) In the case of contact between two beams, formula (A.13) takes the following counterpart form (A.14) Numbering of beams corresponds to a positive sense of a normal to the beam surface. If there is contact, then q k > 0. As it has been already mentioned, the coefficient K is close to one.