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

Nonlinear Dynamics, Nov 2017

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.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

https://link.springer.com/content/pdf/10.1007%2Fs11071-017-3939-6.pdf

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

Contact interaction of two rectangular plates made from different materials with an account of physical nonlinearity J. Awrejcewicz 0 1 2 V. A. Krysko 0 1 2 M. V. Zhigalov 0 1 2 A. V. Krysko 0 1 2 0 A. V. Krysko Cybernetic Institute, National Research Tomsk Polytechnic University , 30 Lenin Avenue, Tomsk , Russia 634050 1 V. A. Krysko 2 J. Awrejcewicz ( 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. Plate; Variation; Iteration; Physical nonlinearity; Contact interaction; Stiffness variable parameters 1 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–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–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. 2 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]: ( 1 ) A1w1(x , y) = q1(x , y) − qk ψ A2w2(x , y) = q2(x , y) − qk ψ where wi (x , y)—vector functions; qi (x , y)—vector of the external load; i = 1, 2—number of a plate coinciding with a positive sense of the normal; Ω1 = {a1(x , y) ≤ z ≤ b1(x , y)}, Ω2 = {a2(x , y) ≤ z ≤ b2(x , y)}—volumes occupied by plates; z = ai (x , y) and z = bi (x , y), (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 = ci (x , y): Ω+ = {ai (x , y) ≤ z ≤ ci (x , y)}, Ω− = i i {ci (x , y) ≤ z ≤ bi (x , y)}. The material in the space Ωi+ is defined by Young’s modulus Ei+(x , y, z, ei+) and Poisson’s coefficient νi+(x , y, z, ei+) corresponding to extension (+), whereas the space Ωi− is defined by Young’s modulus Ei−(x , y, z, ei−) and Poisson’s coefficient νi−(x , y, z, ei−) corresponding to compres ( 2 ) sion (–). Furthermore, ei+/− stands for the deformation intensity. The stress–strain state, owing to reference [29], is defined by non-constant elasticity modulus and Poisson’s coefficient, and the following relations hold σx+/− = σy+/− = τx+y/− = − 1 − νi+/− Ei+/− Ei+/− 1 − νi+/− Ei+/−z 1 + νi+/− 2 εx + νi+/−εy , 2 εy + νi+/−εx , εx y . In physically linear problems, based on the introduced assumptions, the operator Ai (i = 1, 2) takes the following form ∂2wi ∂2wi B11,i ∂ x 2 + B10,i ∂ y2 ∂2 + ∂ y2 ∂2 + 2 ∂ x ∂ y E11,i + (−1)n+1 E10,i E01,i E00,i − E21,i + (−1)n E20,i , Emn,i = bi (x,y) Ei zm 1 + (−1)nνi dz, ai (x,y) n = 0, 1, m = i = 1, 2. The following notation is employed: z = ai (x , y) and/or z = bi (x , y), (x , y) Ωi denote equations of external surfaces of the plates, which allow for inclusion of variable thickness of each plate; Ei (x , y, z, ei(i)), νi (x , y, z, ei(i)) describe different stiffness modulus and Poisson’s coefficients, respectively; ei(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 that Ei (x , y, z, ei(i)),vi (x , y, z, ei(i)) are non-constant parameters dependent on the deformation state in a plate point, and they are governed by the following formulas 9K1i Gi , νi = 21 33KK11ii −+2GGii , Ei = 3K1i + Gi where K1i = const. The theory of deformation yields the following shear modulus where σi(i), ei(i) stand for intensity of stresses and deformations of a plate, respectively, and ek(i) = √32 e(xix) − e(yiy) 2 + e(yiy) − ez(iz) 2 + e(xix) − ez(iz) 2 + 1.5 eix y 2 1/2 . The value of ez(iz) is defined through the condition of a plane stress state (szz = 0): e(xix) + e(yiy) . Two types of boundary conditions are analysed, i.e. νi ez(iz) = − 1 − νi (i) simple support wi ∂Ω = ∂2wi ∂n2 ∂Ω = 0, i = 1, 2; , ( 4 ) Gi = 3 1 σi(i) ei(i) e(i) i , ( 6 ) ( 7 ) ( 8 ) ( 9 ) ( 10 ) ( 11 ) ( 12 ) (ii) rigid clamping wi ∂Ω = ∂wi ∂n ∂Ω = 0, i = 1, 2. Note that in mathematical models ( 1 ), ( 4 ), ( 9 ), ( 10 ) and ( 1 ), ( 5 ), ( 9 ), ( 10 ), the contact pressure is proportional to the transversal crushing w1 − h1 − w2 in a zone of contact between plates, i.e. we have E qk (x , y) = k h (w1 − h1 − w2), and the function ψ has the following form ψ = 1 + sign(w1 − h1 − w2) /2. Here h1 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 h1) 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 A1w1(x , y) = q1(x , y) − k Eh1 (w1 − h1 − w2)ψ, A2w2(x , y) = q2(x , y) + k Eh2 (w1 − h1 − w2)ψ. ( 13 ) 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 ⎧ ⎨⎪⎪⎪⎪⎪⎪ ⎪⎪⎪⎪⎪⎪⎩ A1 w(n) 1 + k Eh ψ(n−1)w1(n) = q1(x , y) + k Eh1 h1 + w(n−1) ψ(n−1), 2 A2 w(n) 2 + k Eh2 w1(n−1) − h1 ψ(n−1). + k Eh ψ(n−1)w2(n) = q2(x , y) ( 14 ) 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 E1, h1/E2, h2, and we have qk (x , y) = 2k E1 h1 E1h2 1 + E2h1 (w1 − h j − w2). ( 15 ) 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. 3 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 Aw (x , y) = q (x y) , (x , y) ∈ Ω = (a, b) × (c, d) , 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). ( 16 ) 3.1 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 ( AwN − q, ϕ j (x ))H(a,b) = 0 j = 1, 2, . . . , N , ( 18 ) Here ( Au N − q, ϕ j (x )) stands for notation of a scalar product in the space H (a, b), H = H (a, b) × H (c, d) is the space in which the exact solution w0 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. 3.2 Method of Vaindiner (MV) N i=1 where ψi (y) and ui (x ) are given functions, ϕi (x ) and vi (x ) are the being searched functions satisfying the system of projected equations of the following form ( 20 ) ( AwN − q, ψ j (y))H(c,d) = 0, ( AwN − q, u j (x ))H(a,b) = 0, j = 1, 2, . . . , N , where H ((a, b) × (c, d)) = H (a, b) × H (c, d) is the space, where the exact solution w0 exists. 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. 3.3 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 ( Au N − q, ϕ j (x ))H(a,b) = 0, ( 21 ) ( Au N − q, ψ j (y))H(c,d) = 0 j = 1, 2, . . . , N , ( 22 ) in the following way. A certain system of N functions with regard to one of the variables, let us say ϕ 0j(x ) j = 1, 2, . . . , N , is given. Then, the system ( 21 ) defines a system of functions ψ (j1)(y), which are substituted into the system ( 22 ), and a new choice with regard to the variable x is defined as ϕ(j2)(x ). Then, using the chosen functions, the system ( 21 ) yields a new function ψ (j3)(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 In this case a solution to Eq. ( 15 ) is assumed to be as follows [34]: wN = ϕi(k−1)(x )ψ (jk)(y), ( 23 ) N i=1 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. 3.4 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 w1 = ϕ1(k−1)(x ) ψ1(k)(y). The new equation is defined as follows Aw2(x , y) = q(x , y) − Aw1(x , y), 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 w2(x , y) = ϕ2(k−1)(x )ψ2(k)(y). The next new equation follows: Aw3(x , y) = q(x , y)− Aw1(x , y)− Aw2(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: w(x , y) = wn(x , y). N n=1 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. ( 24 ) ( 25 ) ( 26 ) ( 27 ) 3.5 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 u0(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 Aw2(x , y) = q(x , y) − Aw1(x , y) ( 29 ) defined through the so far described way, i.e. w2(x , y) = ϕ2(k)(x )ψ2(k−1)(y) + u(2k−1)(x )v(k)(y). If the error between u1 and (u1 +u2) is less than a given quantity, then the process is terminated. If not, then the following equation is constructed Aw3(x , y) = q(x , y) − Aw1(x , y) − Aw2(x , y), ( 30 ) and w3(x , y) is defined, and so on. The final solution to Eq. ( 16 ) is defined by the following series N j=1 wN (x , y) = where ψi (y) and ui (x ) are the given functions, ϕi (x ) and vi (x ) are the searched functions defined by the system of projection Eq. ( 18 ). However, while the found function w1(x , y) is treated as the final one in the MV, it is used to define the operator Aw1(x , y) on the right-hand side of the following equation which is solved by the MV. If the error between w1 and w1 + w2 is smaller than a given quantity, then the iterative process is stopped. Otherwise, we construct the following equations Aw3(x , y) = q(x , y) − Aw1(x , y) − Aw2(x , y), ( 38 ) and we define w3(x , y) with the help of the MKV. As the final solution of Eq. ( 16 ), we take the following 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 w1(x , y) found in the MKV are the final ones, they can be employed to construct the operator Aw1(x , y) on the right-hand side of equation Aw2(x , y) = q(x , y) − Aw1(x , y), ( 33 ) which is solved using the MKV. If the error between w1 and w1 + w2 is smaller than a given quantity, then the iterative process is stopped. If this is not the case, than we construct the following equation Aw3(x , y) = q(x , y) − Aw1(x , y) − Aw2(x , y), ( 34 ) and we define w3(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 ψ 0(y) and u0(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 )(x , y) = ϕ( 1 )(x )ψ ( 1 )(y) + u( 1 )(x )v( 1 )(y). ( 41 ) 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 ( 35 ) ( 39 ) ( 40 ) 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. 3.9 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 Aw(x , y) = D∇2∇2w(x , y) = q(x , y), ( 42 ) which has known exact solution. The following boundary conditions are applied: The space Ω = (0, a) (0, b) has the boundary ∂Ω. The differential operator ( 42 ) and the employed boundary conditions are presented in their counterpart nondimensional form (i) simple support ∂2w w ∂Ω = ∂n2 ∂Ω = 0; (ii) stiff clamping ∂w w ∂Ω = ∂n ∂Ω = 0. ( 43 ) ( 44 ) , 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 ε = wk (x,yw)k−(wx,ky−)1(x,y) < 10−3. The obtained results are reported in Table 1 for q¯ (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. 4 Mathematical background 4.1 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 Aw(x , y) = q(x , y); x , y ∈ Ω(x , y), ( 46 ) where A stands for a certain operator defined on the manifold D( A) of the Hilbert space L2(Ω); 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 ui (x ) and vi (y) are defined by the following system of equations X X ( AwN − q) u1 (x) dx = 0, ( AwN − q) v1 (y) dy = 0, Y ( AwN − q) u N (x) dx = 0, ( AwN − q) vN (y) dy = 0, Y in the following way. A certain system composed of N functions with respect to one of the variables, for instance, u01(x ), u02(x ), . . . , u0N (x ) is given. Then, the first N equations of the system ( 48 ) yield N functions v11(x ), v21(x ), . . . , v1N (x ). Next, the obtained functions are employed to create a new set of functions x − ( 49 ) ( 50 ) u21(x ), u22(x ), . . . , u2N (x ), which is further used to construct a set of new functions with respect to the variable y, i.e. v13(x ), v23(x ), . . . , v3N (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 1 If A is a positively defined operator with its action space D( A) ⊂ HA, then the sequence of elements α = w1k (x , y) − w0 HA is monotonously decreasing, i.e. for arbitrary and j, if only i ≥ j, then w1i − w0 HT ≤ w1j − w0 HT . Theorem 2 Let each element of the basis system of the 0 space W2m (X × Y ) has the form θi (x , y) = ϕi (x )ψi (y), 0 where is a basis system in the space W2m (X × Y ), and, 0 respectively, {Ψi (X × Y )} in the space W2m (Y ), and in order to get an arbitrary -th approximation regarding MVI, the components of the elements of the basis system {θi (X × Y )} are taken as initial conditions. Then, for sufficiently large N , the MVI gives a unique approximate solution, and the sequence {WN } is convergent 0 with regard to the norm of the space W2m (X × Y ), and tends to the exact solution ω0 independently of a number of steps k, which can be defined for each of the N -th approximation, i.e. wkN − w0 W02m , N → ∞ The ordinary differential equations obtained on each step of the iterational procedure by the method of variational iterations can be solved with numerous methods such as reduction of the problem to the algebraic equations by FDM with the successive solutions of the obtained algebraic equations by the Gauss method. 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: w1(x ) = 1; w2(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%. 4.2 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 R2 be an Euclidean plane with a Descartes basis: Ωi R2—area on a plane with the boundary ∂Ωi (i = 1, 2) , Ω¯ i = Ωi ∪ ∂Ωi (x , y) Ωi ; Ω∗—subarea of Ωi , ∀i, Ω∗ ⊆ ΩI ; ni —external normal to ∂Ωi . We consider a space plate-type object composed of two contacting plates (Fig. 2) obeying Kirchhoff’s hypotheses and governed by the following PDEs: A1 (w1(x , y)) + k w1ψ (x , y) = q1 + k Eh1 (w2 + h1)ψ (x , y), A2 (w2(x , y)) + k w2ψ (x , y) E1 h E2 h = q2 (1 − ψ (x , y)) + k Eh2 (w1 − h1)ψ (x , y), with the boundary conditions ( 10 ). The function defining the contact zone Ω∗ of the plates is as follows ψ (x , y) = 1, (x , y) ∈ Ω∗, 0, (x , y) ∈/ Ω∗, and q1(x , y), q2(x , y) are the functions of external load acting on the first and the second plate, respectively; ( 51 ) ( 52 ) the operators Ai (wi ) have the form ( 5 ), which in the case of elastic problems is governed by ( 4 ); · A— norm in a normalized space A; (·, ·)—scalar product in the Hilbert space L2. The employed notation of the functional spaces coincides with that given in reference [40], whereas E1, E2 stand for different moduli in tension and compression. In order to solve the problem ( 51 ), ( 10 ), the following iterational algorithm is employed: A1 w1(n+1)(x , y) + k Eh1 w1(n+1)ψ (x , y) = q1 E1 + k h w(n) 2 + h1 ψ (x , y), A2 w2(n+1)(x , y) + k Eh2 w2(n)ψ (x , y) = q2 (1 − ψ (x , y)) E2 + k h w(n) 1 − h1 ψ (x , y), w(n+1) 1 |∂Ω1 = w(n+1) 2 |∂Ω2 = ∂w1(n+1) Fig. 2 Two plates with a clearance constructions made from materials of different moduli, i.e. materials having different values of the elastic modulus and the Poisson’s coefficients regarding tension and compression. Remark 2 The given theorem implies the theorem reported in reference [39] valid for a material with one modulus. 5 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 σi(i)+(ei(i)+) and σi(i)−(ei(i)−) 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 σv(i) = 13 σ X(i) + σY(i) + σZ(i) , one may construct a diagram σi(i)(ei(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. Fig. 3 The diagram σi (ei ) 5.1 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 σi=σs 1 − exp ei − eis , where ei S = 3σσS0 is a conventional intensity of plasticity deformation, and hence, eis σ = σ1 ei 1 − exp ei − eis . While approximating the σi (ei ) dependence shown in Fig. 3, one may use the following formula σi = σs 1 − exp ei − eis ei 1 + eis . (57) The coefficient b can be identified assuming that the curve L goes through the point (σi(∗), ei(∗)) ∈ L , and hence, ei b = − eis 1 + eeii ∗S ln σi∗ 1 − σs . This allows for more feasible tracing of the velocity of approximating σs by σi∗ with an increase in ei . In the case of a material with different moduli (Fig. 4) and with an account of σi (ei ) in the case of physical nonlinearity, we have (55) (56) (58) The conventionally recognized notation ± for σi , ei , eis 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 × 105 Pa, G1 = G2 = 0.383 × 1011 Pa (here stresses are measured with respect to G1(h/a)2), the non-dimensional deformations e = ei (a/ h)2, whereas the plasticity deformations are described by the formula eS = ei S (a/ h)2. In reference [41], it has been shown that the properties of aluminium are well approximated by the dependencies σi (ei ) through formulas (55), (56) and they are validated experimentally. Owing to the data reported in this reference, we take k = k1 = 2G1, and the σi(i)(ei(i)) function is governed by the following formula ei(i)± − es(i)± σi(i)± = α±σs(i) 1 − exp , (60) the algorithms developed by us for estimating the dependence σi(i)(ei(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 Bi 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, Ei , ν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(w1(0.5, 0.5)) is reported in Fig. 5, q(x,y) a4 where q¯ (x , y) = 12(1−ν12) E1h4 . 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(w1(0.5, 0.5)) are reported in Fig. 6. The value of deflection, for which a contact occurred, is denoted by a dot (w1 = 0.01). In the case of contact problems, the differences exhibited by the curves q(w1) 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 (α− = 1) materials is the same, there is an essential quantitative difference. The surface M¯ x (x , y) = − ∂∂x D1 ∂∂wx 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 surfaces M¯ 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. 6 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 (resonance phenomenon). The wavelet surface w(s, t0) = w(s, t0 exp[ j ϕs (t0)]) characterizes the system behaviour in each time interval s at an arbitrary time instant t0. The quantity w(s, t0) describes the occurrence and intensity of the corresponding time interval s at the time instant t0. The integral energy distribution of the wavelet spectrum with respect to time intervals E (s) = w(s, t0) 2dt0 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) ⎣ ⎣ + + ⎡ E + ci (x,y) z2dz 1 ai (x,y) In a particular case of one-modulus (E1+/− = E2+/− = E and ν1+/− = ν2+/− = ν) physically linear material (E = const, ν = const), PDEs (61) can be reduced to the counterpart non-dimensional form as follows Fig. 7 Surfaces Mx (x, y) for different (a) and the same (b) modulus corresponding to the upper and lower plate, respectively Fig. 8 Surfaces Mx (x, y) for the upper plate with α− = 0.8 (a) and α− = 1.0 (b), respectively ⎧⎪ 12(11−ν2) ∇λ4w1 + ∂∂2tw21 + ε ∂∂wt1 ⎪⎪⎨ −K (w1 − w2 − h∗)Ψ = q(t ), ⎪⎪ 12(11−ν2) ∇λ4w2 + ∂∂2tw22 + ε ∂∂wt2 ⎩⎪ +K (w1 − w2 − h∗)Ψ = 0. System of Eq. (62) has been obtained by the following introduced relations between dimensional and non-dimensional parameters: x = ax , y = a y¯; q = q Ea(2hb)24 , τ = ahb " Eγg , λ = ab , where a, b—plate dimension with respect to x and y, respectively; t — time; ε—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: wm = 0; wm |x = 0; for x = 0; 1, wm = 0; wm |y = 0; for y = 0; 1; m = 1, 2. The initial conditions are as follows wm (x , y)|t=0 = 0, w˙ m |t=0 = 0, m = 1, 2. We take the control function ψ = 21 [1 + sign(w1 − h ∗ −w2)]; Ψ = 1, if w1 > w2 + hk , i.e. there is contact between plates; otherwise Ψ = 0; w1, w2— (62) (63) (64) 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 ) = q0 sin(ω0t ). 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 wi , 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 wm = N N i=1 j=1 Ai(mj)(t ) sin(i π x ) sin( j π y), m = 1, 2 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: q0 = 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 (65) 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 history 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 wn(x , y), where a solution to the dynamical problem can be presented in the following form w(x , y, t ) = A(t )wn(x , y), (66) 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. 7 Concluding remarks The following general conclusions can be formulated as a result of the carried out research. Acknowledgements This work has been supported by the Ministry of Education and Science of the Russian Federation by the Grant No 2.1642.2017. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Appendix A Coupling of contact pressure and transverse clamping of a thin beam Let us show that in order to take into account the transverse clamping, one can employ the Winkler’s hypothesis. To solve this problem, we begin with a solution to q x Integrating (A.2) with respect to z, and taking into account that 0l1 σx dz = 0 and since there is no load in the x direction, we obtain We take q(x ) = −σz (x , l1). In the case of a plane in a deformable state, we have εz = ∂uz (x , z) ∂ z (q − υ2)σz − υ (1 + υ) σz . 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. (A.1) (A.2) 2l1 2 l Fig. 10 Contact between a thin layer of the thickness 2l1 and a die the problem of contact of a layer of thickness 2l1 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 q(x ) = Am cos ∞ m=0 mπ x l . The stresses are derived in reference [46], and they take the following form Am Fm−(z) cos β x , Am Fm+(z) cos β x , (βl1chβl1 ± shβl1)chβ z − β zshβ zshβl1 sh2βl1 + 2βl1 , σx = 2 σz = 2 Fm± = β = mπ l . 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 uz (x , l1) is given on the interval −l ≤ x ≤ l, satisfies the following equation (A.4) (A.5) q(ξ )K dξ = π s0uz (x , l1), (A.6) ∞ −∞ 2l z Fig. 11 Contact between a thin layer and two symmetrically located dies Solution to Eq. (A.6) in the form was reported in reference [48], where the approximation L(uu) = A1 + #i∞=1 Bi u(2i)( p) was used to obtain the following relation: K ∗( p) = πA δ0( p) + π #i∞=1 (−1)i Bi δ(2i)( p) where δ(2i) stands for the 0 derivative of 2i -th order (Fig. 11). In what follows, we find the solution in an explicit form of the mentioned formulas [48]: = l1δ0(x − ξ ), δ0(x − ξ ) l −l l −l δ0( p) = δ0 x − ξ l1 = −δ0(ξ − x ), f (ξ )δ0(2i)(x − ξ )dξ = f (2i)(x ), |x | ≤ l. (A.9) Owing the evenness of the function K ∗( p), we have K ∗ ξ −l1x = K ∗ xl−1ξ . Substituting the series of u/L(u) into the solution of Eq. (A.6), carrying out the integration in the interval −l, l, and taking into account the already known relations for the generalized functions, we get 1 E q(x ) = 2(1 − υ2) l1 uz (x , l1) − B1l12uz + B2l14uzI V − · · · . (A.10) x (A.8) l 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 3(1−υ) the case of plates yields qk = 4(1+υ)(1−2υ) Eh Δ + o(h2 Δ0w) . In the case of beams, taking into account t43heEh PΔoi+ssoon(’hs2Δco0ewffi)c.ient υ = 0, it yields qk = 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–54], and the following simple formula holds qk = K (w − a) , a > 0. qk = K (w1 − a − w2). h In the case of contact between two beams, formula (A.13) takes the following counterpart form h Numbering of beams corresponds to a positive sense of a normal to the beam surface. If there is contact, then qk > 0. As it has been already mentioned, the coefficient K is close to one. (A.13) (A.14) E E 1. Jones , R.M.: Apparent flexural modulus and strength of multimodulus materials . J. Compos. Mater . 10 , 342 - 354 ( 1976 ) 2. Guo , Z.H. , Zhang , X.Q. : Investigation of complete stressdeformation curves for concrete in tension . ACI Mater. J . 84 ( 4 ), 278 - 285 ( 1987 ) 3. Haimson , B.C. , Tharp , T.M.: Stresses around boreholes in bilinear elastic rock . Soc. Petrol. Eng . 14 ( 2 ), 145 - 151 ( 1974 ) 4. Gilbert , G.N.J.: Stress/strain properties of cast iron and Poisson's ratio in tension and compression . Brit. Cast Res. Assn. J . 9 , 347 - 363 ( 1961 ) 5. Patel , H.P. , Turner , J.L. , Walter , J.D. : Radial tire cordrubber composites . Rubber Chem. Technol . 49 ( 4 ), 1095 - 1110 ( 1976 ) 6. Barak , M.M. , Currey , J.D. , Weiner , S. , Shahar , R.: Are tensile and compressive Young's moduli of compact bone different . J. Mech. Behav. Biomed . 2 ( 1 ), 51 - 60 ( 2009 ) 7. Bertoldi , K. , Bigoni , D. , Drugan , W.J.: Nacre: an orthotropic and bimodular elastic material . Compos. Sci. Technol . 68 , 1363 - 1375 ( 2008 ) 8. Geim , A.K. : Graphene: status and prospects . Science 324 ( 5934 ), 1530 - 1534 ( 2009 ) 9. Tsoukleri , G. , Parthenios , J. , Papagelis , K. , et al.: Subjecting a graphene monolayer to tension and compression . Small 5 ( 21 ), 2397 - 2402 ( 2009 ) 10. Yao , W. , Ma, J. , Gao , J. , Qiu , Y. : Nonlinear large deflection buckling analysis of compression rod with different moduli . Struct. Eng. Mech . 54 ( 5 ), 855 - 875 ( 2015 ) 11. Yao , W.J. , Ye , Z.M. : Analytical method and the numerical model of temperature stress for the structure with different modulus . J. Mech. Strength 27 ( 6 ), 808 - 814 ( 2005 ) 12. Bert , C.W. : Models for fibrous composites with different properties in tension and compression . ASME J. Eng. Mater. Technol . 99 , 344 - 349 ( 1977 ) 13. Reddy , J.N. , Chao , W.C. : Finite-element analysis of laminated bimodulus composite-material plates . Comput. Struct . 12 , 245 - 251 ( 1980 ) 14. Bruno , D. , Lato , S. , Zinno , R.: Nonlinear analysis of doubly curved composite shells of bimodular material . Compos. B 3 , 419 - 435 ( 1993 ) 15. Zinno , R. , Greco , F. : Damage evolution in bimodular laminated composites under cyclic loading . Compos. Struct . 53 , 381 - 402 ( 2001 ) 16. Patel , B.P. , Lele , A.V. , Ganapathi , M. , Gupta , S.S. , Sambandam , C.T.: Thermo-flexural analysis of thick laminates of bimodulus composite materials . Compos. Struct . 63 , 11 - 20 ( 2004 ) 17. Ambartsumyan , S.A. , Wu , R.F. , Zhang , Y.Z. : Elasticity Theory of Different Moduli. China Railway Publishing House , Beijing ( 1986 ) 18. Jones , R.M. : Stress-strain relations for materials with different moduli in tension and compression . AIAA J . 15 ( 1 ), 16 - 23 ( 1977 ) 19. Zhang, Liang, Zhang, Huiting, Jian, Wu, Yan, Bo, Mengkai, Lu: Parametric variational principle for bi-modulus materials and its application to nacreous bio-composites . Int. J. Appl. Mech . 8 ( 6 ), 1650082 ( 2016 ) 20. Hertz , H.: Gesammelte Werke. Johann Ambrosius Barth , Leipzig ( 1895 ) 21. Duvaut , G. , Lions , J.-L.: Les in Equations en Mecanique et en Physique . Dunod, Paris ( 1972 ) 22. Banks , H.T. , Hu , S. , Kenz , Z.R.: A brief review of elasticity and viscoelasticity for solids . Adv. Appl. Math. Mech. 3 , 1 - 51 ( 2011 ) 23. Migorski , S. , Ochal , A. , Sofonea , M. : Nonlinear Inclusions and Hemivariational Inequalities . Models and Analysis of Contact Problems . Springer, Berlin ( 2013 ) 24. Sofonea , M. , Matei , A. : Mathematical Models in Contact Mechanics . Cambridge University Press, Cambridge ( 2012 ) 25. Migorski , M. , Ochal , A. , Sofonea , M.: Analysis of a quasistatic contact problem for piezoelectric materials . J. Math. Anal. Appl . 382 , 701 - 713 ( 2011 ) 26. Migorski , S. , Ochal , A. , Sofonea , M.: A dynamic frictional contact problem for piezoelectric materials . J. Math. Anal. Appl . 361 , 161 - 176 ( 2010 ) 27. Kantor , B.Y. : Nonlinear Problems of the Theory of Nonhomogeneous Shallow Shells . Naukova Dumka , Kiev ( 1971 ). in Russian 28. Zhao , H.L. , Ye , Z.M.: Analytic elasticity solution of bimodulus beams under combined loads . Appl. Math. Mech . 36 , 427 - 438 ( 2015 ) 29. He , X.T. , Chen , Q. , Sun , J.Y. , Zheng , Z.L. , Chen , S.L. : Application of the Kirchhoff hypothesis to bending thin plates with different moduli in tension and compression . J. Mech. Mater. Struct . 5 , 755 - 769 ( 2010 ) 30. Birger , I.A. , Mavlyutov , R.R. : The Resistance of Materials. Nauka , Moscow ( 1986 ). in Russian 31. Kantorovich , L.V. : One direct method of approximate solution of the problem of minimum of the double integral . News Acad. Sci. SSSR 5 , 647 - 652 ( 1933 ). in Russian 32. Vlasov , V.Z.: A new practical method for calculating folded coatings and shells . Constr. Ind . 11 , 33 - 37 ( 1932 ) 33. Fletcher , C.A.J. : Computational Galerkin Methods . Springer, Berlin ( 1984 ) 34. Vaindiner , A.I. : On a generalization of the BubnovGalerkin-Kantorovich method of approximate solution of boundary value problems . Bull. Moscow St. Univ. Math. Mech. 2 , 1001 - 1003 ( 1967 ). in Russian 35. Kirichenko , V.F. , Krysko , V.A. : Substantiation of the variational iteration method in the theory of plates . Sov. Appl. Mech . 17 ( 4 ), 366 - 370 ( 1981 ). in Russian 36. Arganovskyy , M.L. , Baglay , R.D.: On an expansion in Hilbert space and its applications . Comput. Math. Math. Phys . 17 ( 4 ), 871 - 878 ( 1977 ). in Russian 37. Baglay , R.D., Smirnov , K.K. : To the processing of twodimensional signals on a computer . Comput. Math. Math. Phys . 15 ( 1 ), 241 - 248 ( 1975 ). in Russian 38. Mikhlin , S.G. : Variational Methods in Mathematical Physics. International Series of Monographs in Pure and Applied Mathematics , vol. 50 , p. 584 . Pergamon Press, Oxford ( 1964 ) 39. Krysko , A.V. , Awrejcewicz , J. , Zhigalov , M.V. , Krysko , V.A. : On the contact interaction between two rectangular plates . Nonlinear Dyn . 84 ( 4 ), 2729 - 2748 ( 2016 ) 40. Ladyzhenskaya , O.A. , Uraltseva , N.N. : Linear and Quasilinear Elliptic Equations . Academic Press, New York ( 1968 ) 41. Ohashi , Y. , Murakami , S.: The elasto-plastic bending of a clamped thin circular plate . In: Proceedings of the Eleventh International Congress of Applied Mechanics Munich , pp. 212 - 223 ( 1964 ) 42. Chen , X. , Xiang , J. , Li , B. , He , Z. : A study of multiscale wavelet-based elements for adaptive finite elements analysis . Adv. Eng. Soft . 41 , 196 - 205 ( 2010 ) 43. Lepik , U. : Solving integral and differential equations by the aid of nonuniform Haar wavelets . Appl. Math. Comput . 198 , 326 - 332 ( 2008 ) 44. Awrejcewicz , J. , Krysko , V.A. , Papkova , I.V. , Krysko , A.V. : Deterministic Chaos in One Dimensional Continuous Systems . World Scientific, Singapore ( 2016 ) 45. Awrejcewicz , J. , Krysko , V.A. : Nonclassical Thermoelastic Problems in Nonlinear Dynamics of Shells. Springer, Berlin( 2003 ) 46. Vorovich , I.I. , Aleksandrov , V.M. , Babenko , V.A. : Nonclassical Hybrid Problems of Theory of Elasticity . Nauka, Moscow ( 1974 ). in Russian 47. Aleksandrov , V.M. , Mkhitarian , S.M. : Contact Problems of Bodies with Thin Covers and Layers . Nauka, Moscow ( 1983 ). in Russian 48. Korn , G. , Korn , T. : Guide for Mathematics . Nauka, Moscow ( 1968 ) 49. Bogatyrenko , T.D., Kantot , B.Y. , Lipovskiy , D.E. : On the contact of rotational shells via thin elastic layer . Bull. Inst. Mach. Des. AN SSSR 8 , 11 - 49 ( 1982 ). (in Russian) 50. Bloch , V.M. : On the model choice in the contact problems of thin-walled bodies . Appl. Mech . 13 , 34 - 42 ( 1977 ). (in Russian) 51. Hencky , H.: Zur Theorie plastischer Deformationen und der hierdurch im Material hervorgerufenen Nachspannungen . ZAMM 4 , 323 - 334 ( 1924 ) 52. Galin , L.A. : Development of Theory of Contact Problems in SSSR . Nauka, Moscow ( 1976 ). (in Russian) 53. Grigoluk , E.I. , Tolmachev , V.M.: Cylindrical bending of a plated by rigid stamps . Appl. Math. Mech . 39 , 876 - 883 ( 1975 ). (in Russian) 54. Bloch , M.V. : One-dimensional one sided contact of rods, plates, and shells. Theory of shells and plates . In: Proceedings of the 9th Soviet Conference on Theory of Plates and Shells. Leningrad , pp. 25 - 28 ( 1975 ) (in Russian)


This is a preview of a remote PDF: https://link.springer.com/content/pdf/10.1007%2Fs11071-017-3939-6.pdf

J. Awrejcewicz, V. A. Krysko, M. V. Zhigalov, A. V. Krysko. Contact interaction of two rectangular plates made from different materials with an account of physical nonlinearity, Nonlinear Dynamics, 2017, 1191-1211, DOI: 10.1007/s11071-017-3939-6