Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling hp-Version of Collocation and Least Residuals Method in Mechanics of Laminated Composite Plates Sergey Golushko1,2 , Semyon Idimeshev2 , and Vasiliy Shapeev3 1 Institute of Computational Technologies Siberian Branch of the Russian Academy of Sciences, Academician M.A. Lavrentjev ave. 6, Novosibirsk, Russia 2 Design Technological Institute of Digital Techniques Siberian Branch of the Russian Academy of Sciences, Academician Rzhanova ave. 6, Novosibirsk, Russia 3 Institute of Theoretical and Applied Mechanics Siberian Branch of the Russian Academy of Sciences, Institutskaya Str. 4/1, Novosibirsk, Russia {s.k.golushko,idimeshev,vshapeev}@gmail.com http://www.ict.nsc.ru/en Abstract. A version of collocations and least residuals method (CLS) based on polynomial approximation of high degree (p - approach) was proposed and implemented. In rectangular domains collocation points are selected using the roots of Chebyshev polynomials and approximate so- lution is represented in the form of direct products of Chebyshev polyno- mials series. It was shown that the use of p - approach in the CLS method allows to obtain numerical solutions with high accuracy and to implement complex boundary conditions with no special techniques. The numeri- cal method used to solve a problem of bending of laminated anisotropic rectangular plates within frameworks of classical laminated plate the- ory, first order shear deformation theory and Grigolyuk-Chulkov’s broken line theory. Several specific example problems are solved, including fixed three-ply laminates with transversely isotropic layers under transverse uniform loading. Keywords: Collocations and least residuals method, Chebyshev poly- nomials, plate theory, spectral methods, composite materials 1 Introduction The collocation and least residuals method (CLS) is an efficient method for numerial solution of boundary value problems both for systems of ordinary and partial differential equations. It is based on the collocation method (CM) [1], with approximate solution is represented as a linear combination of basis functions in some functional space. To determine it unknown coefficients in CM residual of equations 𝑅(π‘₯) vanishes at given points (collocation points) 𝑅(π‘₯π‘π‘œπ‘™ 𝑖 ) = 0, {π‘₯π‘π‘œπ‘™ } – collocation points. (1) 299 Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling The main difference between CLS method and CM is the minimizing technique of 𝑅(π‘₯). In CLS method we minimize some functional of residual in the collocation points [2, 3], instead of the condition (1). The CLS method is used to minimize residual in 𝐿2 norm βˆ‘οΈ ‖𝑅(π‘₯π‘π‘œπ‘™ 2 𝑖 )β€–2 β†’ min . (2) 𝑖 In CLS method the number of equations can exceed the number of unknown coefficients in representation of the solution. The solutions of arising overde- termined systems of linear algebraic equations (SLAE) are defined in the sense of (2) (least squares). In comparison with CM (1) the obtained overdetermined SLAE is often better conditioned and leads to less nonphysical oscillations in numerical solutions. Similar regularization approaches are applied in the finite element method (Least squares finite element methods) [4]. In this paper, an approximate solution is represented as a linear combination of polynomials of high degrees (p - approach) that is typical for spectral methods. This allows to obtain numerical solutions of high accuracy at low computational cost. Term (2) makes the implementation of p - approach more convenient in CLS method when compared with CM. This modification based on p - approach is called hp - version of CLS method. We will demonstrate the application of a method to solving problems of solid mechanics – bending of laminated anisotropic rectangular plates. On a practical level, for calculating the stress and displacement fields of such structures the the- ories of plates are used. They lead to a smaller computational efforts compared to three-dimentional elasticity formulation. Boundary value problems arising in a plate theories have a number of features that present difficulties for many well-known numerical methods. First, govern- ing equations of plate theories may contain derivatives of high orders. Second, boundary conditions may be quite complicated, for example, in a form of linear combination of functions and their higher order derivatives. Third, equations of the plate theory may contain small parameters in the derivatives. These features cause serious difficulties for widely used finite differences and finite element meth- ods. The use of both p - approach and term (2) in hp - version of CLS method may resolve these difficulties and obtain high accuracy solutions at relatively low computational efforts. 2 Formulation of Problem and Governing Equations Let us consider a static bending of laminate composed of 3 layers of constant thickness (Fig. 1). Layers are transversely isotropic with material symmetry axis in the plate’s plane. Layers orientation scheme is πœƒ1 = πœƒ3 = 0, πœƒ2 = πœ‹/2, where πœƒπ‘˜ is an angle measured counterclockwise from the π‘₯ coordinate axis to the π‘˜-th layer material symmetry axis. 300 Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling z z q(x,y) 3 h z3 2 x h z2 x y 2 z1 1 b a Fig. 1. Rectangular multilayered plate under transverse loading; π‘Ž, 𝑏, β„Ž β€” plate’s di- mensions in the directions π‘₯, 𝑦, 𝑧 respectively, π‘§π‘˜ β€” π‘˜-th layer lower surface coordinate, π‘˜ = 1, 2, 3. Engineering constants of transversely isotropic material are [5] 𝐸𝐿 = 25 Mpsi, 𝐸𝑇 = 1 Mpsi, 𝐺𝐿𝑇 = 0.5 Mpsi, (3) 𝐺𝑇 𝑇 = 0.2 Mpsi, πœˆπΏπ‘‡ = πœˆπ‘‡ 𝑇 = 0.25. Here 𝐸, 𝐺, 𝜈 are elasticity and shear modulus, Poisson ratios. 𝐿 signifies the material symmetry axe, 𝑇 the transverse direction. Layers thicknesses are β„Ž1 = β„Ž3 = β„Ž/4, β„Ž2 = β„Ž/2. The upper surface of the plate is under uniform transverse load π‘ž0 , the lower surface is free, and a continuity condition of displacements 𝑒, 𝑣, 𝑀 and stresses πœŽπ‘§π‘§ , 𝜎π‘₯𝑧 , πœŽπ‘¦π‘§ is used on interface surfaces. The corresponding boundary condi- tions are defined on the boundary of the plate. The task is to calculate the stress and displacement fields of such plates. Calculation of thin laminated anisotropic structures within framework of the three-dimensional elasticity is associated with high computational efforts. There- fore, many researchers frequently make use of more robust plate theories, that allows to reduce the dimension of the original problem by excluding the direction of the coordinate 𝑧 from consideration. We consider three plate theories, where transverse shear stresses are simulated differently: classical laminated plate the- ory (CLPT), first order shear deformation theory (FSDT) or Timoshenko’s plate theory [6] and Grigolyuk-Chulkov’s broken line theory (GCT) [7]. Classical laminated plate theory uses the classical Kirchhoff assumption which implies the geometric relationships in the form of πœ•π‘’0 πœ• 2 𝑀0 πœ•π‘£0 πœ• 2 𝑀0 𝑒π‘₯π‘₯ = βˆ’π‘§ , 𝑒 𝑦𝑦 = βˆ’ 𝑧 , 𝑒𝑧𝑧 = 0, πœ•π‘₯ πœ•π‘₯2 πœ•π‘¦ πœ•π‘¦ 2 (οΈ‚ )οΈ‚ 1 πœ•π‘’0 πœ•π‘£0 πœ• 2 𝑀0 𝑒π‘₯𝑦 = + βˆ’π‘§ , 𝑒π‘₯𝑧 = 0, 𝑒𝑦𝑧 = 0. 2 πœ•π‘¦ πœ•π‘₯ πœ•π‘₯πœ•π‘¦ 301 Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling Here 𝑒𝑖𝑗 – strains; 𝑒0 (π‘₯, 𝑦), 𝑣0 (π‘₯, 𝑦), 𝑀0 (π‘₯, 𝑦) – central plane displacements. Constitutive equations for π‘˜-th layer are expressed by βŽ› π‘˜ ⎞ βŽ› π‘˜ π‘˜ π‘˜ βŽžβŽ› ⎞ 𝜎π‘₯π‘₯ 𝑄11 𝑄12 𝑄16 𝑒π‘₯π‘₯ ⎜ π‘˜ ⎟ ⎜ π‘˜ π‘˜ π‘˜ ⎟⎜ ⎟ βŽπœŽπ‘¦π‘¦ ⎠ = βŽπ‘„12 𝑄22 𝑄26 ⎠ βŽπ‘’π‘¦π‘¦ ⎠ , (4) π‘˜ 𝜎π‘₯𝑦 π‘„π‘˜16 π‘„π‘˜26 π‘„π‘˜66 𝑒π‘₯𝑦 βŽ› π‘˜ π‘˜ π‘˜ ⎞ 𝑄11 𝑄12 𝑄16 ⎜ π‘˜ π‘˜ π‘˜ ⎟ π‘˜ π‘˜ π‘˜ 𝑇 βŽπ‘„12 𝑄22 𝑄26 ⎠ = 𝐷1 𝐢1 (𝐷1 ) , π‘„π‘˜16 π‘„π‘˜26 π‘„π‘˜66 where βŽ› π‘˜ π‘˜ ⎞ βŽ› ⎞ 𝐢11 𝐢12 0 cos2 πœƒπ‘˜ sin2 πœƒπ‘˜ βˆ’ sin 2πœƒπ‘˜ ⎜ π‘˜ π‘˜ ⎟ 𝐢1π‘˜ = ⎝𝐢12 𝐢22 0 ⎠ , 𝐷1π‘˜ = ⎝ sin2 πœƒπ‘˜ cos2 πœƒπ‘˜ sin 2πœƒπ‘˜ ⎠ 0 0 𝐢66 π‘˜ (sin 2πœƒ )/2 βˆ’(sin 2πœƒ )/2 cos 2πœƒπ‘˜ π‘˜ π‘˜ π‘˜ Coefficients 𝐢𝑖𝑗 express in terms of the engineering constant as follows: π‘˜ πΈπΏπ‘˜ π‘˜ πΈπ‘‡π‘˜ 𝐢11 = π‘˜ πœˆπ‘˜ , 𝐢22 = π‘˜ πœˆπ‘˜ , 1 βˆ’ πœˆπΏπ‘‡ 𝑇𝐿 1 βˆ’ πœˆπΏπ‘‡ 𝑇𝐿 π‘˜ π‘˜ πœˆπΏπ‘‡ πΈπ‘‡π‘˜ π‘˜ 𝐢12 = π‘˜ πœˆπ‘˜ , 𝐢66 = 𝐺𝐿𝑇 . 1 βˆ’ πœˆπΏπ‘‡ 𝑇𝐿 It is convenient to define the following quantities 3 ∫︁ π‘§π‘˜+1 βˆ‘οΈ 3 βˆ‘οΈ 𝐴𝑖𝑗 = π‘„π‘˜π‘–π‘— 𝑑𝑧 = π‘„π‘˜π‘–π‘— (π‘§π‘˜+1 βˆ’ π‘§π‘˜ ), π‘˜=1 π‘§π‘˜ π‘˜=1 (5) 3 ∫︁ π‘§π‘˜+1 βˆ‘οΈ 3 βˆ‘οΈ 1 𝐷𝑖𝑗 = π‘„π‘˜π‘–π‘— 𝑧 2 𝑑𝑧 = π‘„π‘˜π‘–π‘— (π‘§π‘˜+1 3 βˆ’ π‘§π‘˜3 ), π‘§π‘˜ 3 π‘˜=1 π‘˜=1 where π‘§π‘˜ – coordinates of layers lower surface (Fig. 1). Finaly governing equations for considered problem within CLPT framework are given by πœ• 2 𝑣0 πœ• 2 𝑒0 πœ• 2 𝑒0 (𝐴12 + 𝐴66 ) + 𝐴11 2 + 𝐴66 = 0, πœ•π‘₯ πœ•π‘¦ πœ•π‘₯ πœ•π‘¦ 2 πœ• 2 𝑒0 πœ• 2 𝑣0 πœ• 2 𝑣0 (𝐴12 + 𝐴66 ) + 𝐴22 2 + 𝐴66 = 0, (6) πœ•π‘₯ πœ•π‘¦ πœ•π‘¦ πœ•π‘₯2 πœ• 4 𝑀0 πœ• 4 𝑀0 πœ• 4 𝑀0 (2𝐷12 + 4𝐷66 ) 2 2 + 𝐷11 4 + 𝐷22 = π‘ž0 . πœ•π‘₯ πœ•π‘¦ πœ•π‘₯ πœ•π‘¦ 4 We consider two kinds of boundary conditions: 302 Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling – clamped πœ•π‘€0 π‘₯=0: 𝑒0 = 0, 𝑣0 = 0, 𝑀0 = 0, = 0; πœ•π‘₯ πœ•π‘€0 π‘₯=π‘Ž: 𝑒0 = 0, 𝑣0 = 0, 𝑀0 = 0, = 0; πœ•π‘₯ πœ•π‘€0 𝑦=0: 𝑒0 = 0, 𝑣0 = 0, 𝑀0 = 0, = 0; πœ•π‘¦ πœ•π‘€0 𝑦=π‘Ž: 𝑒0 = 0, 𝑣0 = 0, 𝑀0 = 0, = 0; πœ•π‘¦ – simply-supported πœ•π‘’0 πœ• 2 𝑀0 π‘₯=0: = 0, 𝑣0 = 0, 𝑀0 = 0, = 0; πœ•π‘₯ πœ•π‘₯2 πœ•π‘’0 πœ• 2 𝑀0 π‘₯=π‘Ž: = 0, 𝑣0 = 0, 𝑀0 = 0, = 0; πœ•π‘₯ πœ•π‘₯2 2 πœ•π‘£0 πœ• 𝑀0 𝑦=0: 𝑒0 = 0, = 0, 𝑀0 = 0, = 0; πœ•π‘¦ πœ•π‘¦ 2 πœ•π‘£0 πœ• 2 𝑀0 𝑦=π‘Ž: 𝑒0 = 0, = 0, 𝑀0 = 0, = 0. πœ•π‘¦ πœ•π‘¦ 2 First order shear deformation theory allows transverse shear in a first approximation by defining independent function of rotation of the transverse normal about central surface: πœ‘π‘₯ (π‘₯, 𝑦) and πœ‘π‘¦ (π‘₯, 𝑦). The strains are obtained by πœ•π‘’0 πœ•πœ‘π‘₯ πœ•π‘£0 πœ•πœ‘π‘¦ 𝑒π‘₯π‘₯ = +𝑧 , 𝑒𝑦𝑦 = +𝑧 , πœ•π‘₯ πœ•π‘₯ πœ•π‘¦ πœ•π‘¦ (οΈ‚ )οΈ‚ (οΈ‚ )οΈ‚ πœ•π‘’0 πœ•π‘£0 πœ•πœ‘π‘₯ πœ•πœ‘π‘¦ 𝑒π‘₯𝑦 = + +𝑧 + , πœ•π‘¦ πœ•π‘₯ πœ•π‘¦ πœ•π‘₯ πœ•π‘€0 πœ•π‘€0 𝑒π‘₯𝑧 = + πœ‘π‘₯ , 𝑒𝑦𝑧 = + πœ‘π‘¦ , 𝑒𝑧𝑧 = 0. πœ•π‘₯ πœ•π‘¦ Constitutive equations of FSDT are obtained by adding to (4) expressions for the shear stresses (οΈƒ )οΈƒ (οΈƒ )οΈƒ (οΈƒ )οΈƒ π‘˜ πœŽπ‘¦π‘§ π‘„π‘˜44 π‘„π‘˜45 𝑒𝑦𝑧 π‘˜ = , 𝜎π‘₯𝑧 π‘„π‘˜45 π‘„π‘˜55 𝑒π‘₯𝑧 (οΈƒ )οΈƒ (οΈ‚ )οΈ‚ (οΈ‚ )οΈ‚ π‘„π‘˜44 π‘„π‘˜45 π‘˜ 𝐢44 0 cos πœƒπ‘˜ sin πœƒπ‘˜ = 𝐷2π‘˜ 𝐢2π‘˜ (𝐷2π‘˜ )𝑇 , 𝐢2π‘˜ = π‘˜ , 𝐷2π‘˜ = , π‘„π‘˜45 π‘„π‘˜55 0 𝐢55 βˆ’ sin πœƒπ‘˜ cos πœƒπ‘˜ where stiffness coefficients are express by engineering constants π‘˜ π‘˜ 𝐢44 = 𝐺𝑇 𝑇 , 𝐢55 = 𝐺𝐿𝑇 . 303 Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling Governing equations for FSDT can be written as πœ• 2 𝑣0 πœ• 2 𝑒0 πœ• 2 𝑒0 (𝐴12 + 𝐴66 ) + 𝐴11 + 𝐴 66 = 0, πœ•π‘₯ πœ•π‘¦ πœ•π‘₯2 πœ•π‘¦ 2 πœ• 2 𝑒0 πœ• 2 𝑣0 πœ• 2 𝑣0 (𝐴12 + 𝐴66 ) + 𝐴22 2 + 𝐴66 = 0, πœ•π‘₯ πœ•π‘¦ πœ•π‘¦ πœ•π‘₯2 πœ•πœ‘π‘¦ πœ• 2 𝑀0 πœ•πœ‘π‘₯ πœ• 2 𝑀0 βˆ’π΄44 βˆ’ 𝐴44 βˆ’ 𝐴 55 βˆ’ 𝐴 55 = π‘ž0 , πœ•π‘¦ πœ•π‘¦ 2 πœ•π‘₯ πœ•π‘₯2 πœ• 2 πœ‘π‘¦ πœ•π‘€0 πœ• 2 πœ‘π‘₯ πœ• 2 πœ‘π‘₯ (𝐷12 + 𝐷66 ) βˆ’ 𝐴55 πœ‘π‘₯ βˆ’ 𝐴55 + 𝐷11 + 𝐷66 = 0, πœ•π‘₯ πœ•π‘¦ πœ•π‘₯ πœ•π‘₯2 πœ•π‘¦ 2 πœ• 2 πœ‘π‘₯ πœ•π‘€0 πœ• 2 πœ‘π‘¦ πœ• 2 πœ‘π‘¦ (𝐷12 + 𝐷66 ) βˆ’ 𝐴44 πœ‘π‘¦ βˆ’ 𝐴44 + 𝐷22 2 + 𝐷66 = 0, πœ•π‘₯ πœ•π‘¦ πœ•π‘¦ πœ•π‘¦ πœ•π‘₯2 where coefficients are defined by (5). Boundary conditions for clamped edges in FSDT are written as follows: 𝑒0 = 0, 𝑣0 = 0, 𝑀0 = 0, πœ‘π‘₯ = 0, πœ‘π‘¦ = 0, (π‘₯, 𝑦) ∈ πœ•π›Ί. More details of CLPT and FSDT plate theories are described in [6]. Grigolyuk-Chulkov’s theory is layerwise theory, where mechanical prop- erties of each layer are considered separately. For this purpose, in each layer the rotations of transverse normal about central surface πœ‘π‘˜π‘₯ (π‘₯, 𝑦) and πœ‘π‘˜π‘¦ (π‘₯, 𝑦) are defined. It can be assumed that the GCT is a generalization of the FSDT that takes into account transverse shear stresses in each layer separately. Expressions for geometrical equations in GCT have the form 3 πœ•π‘’0 βˆ‘οΈ πœ•πœ‘π‘–π‘₯ πœ•πœ‘π‘˜ π‘’π‘˜π‘₯π‘₯ = + π‘ƒπ‘˜π‘– + (𝑧 βˆ’ π‘§π‘˜βˆ’1 ) π‘₯ , πœ•π‘₯ 𝑖=1 πœ•π‘₯ πœ•π‘₯ 3 πœ•π‘£0 βˆ‘οΈ πœ•πœ‘π‘–π‘¦ πœ•πœ‘π‘˜π‘¦ π‘’π‘˜π‘¦π‘¦ = + π‘ƒπ‘˜π‘– + (𝑧 βˆ’ π‘§π‘˜βˆ’1 ) , πœ•π‘¦ 𝑖=1 πœ•π‘¦ πœ•π‘¦ (οΈ‚ )οΈ‚ βˆ‘οΈ3 (οΈƒ )οΈƒ (οΈƒ )οΈƒ π‘˜ πœ•π‘£0 πœ•π‘’0 πœ•πœ‘π‘–π‘¦ πœ•πœ‘π‘–π‘₯ πœ•πœ‘π‘˜π‘¦ πœ•πœ‘π‘˜π‘₯ 𝑒π‘₯𝑦 = + + π‘ƒπ‘˜π‘– + + (𝑧 βˆ’ π‘§π‘˜βˆ’1 ) + , πœ•π‘₯ πœ•π‘¦ 𝑖=1 πœ•π‘₯ πœ•π‘¦ πœ•π‘₯ πœ•π‘¦ πœ•π‘€ πœ•π‘€ π‘’π‘˜π‘₯𝑧 = πœ‘π‘˜π‘₯ + , π‘’π‘˜π‘¦π‘§ = πœ‘π‘˜π‘¦ + , π‘’π‘˜π‘§π‘§ = 0, πœ•π‘₯ πœ•π‘¦ where βŽ› ⎞ 0 0 0 𝑃 = βŽβ„Ž1 0 0⎠ . β„Ž1 β„Ž2 0 Boundary conditions for clamped edges in GCT are (π‘˜ = 1, 2, 3) 𝑒0 = 0, 𝑣0 = 0, 𝑀0 = 0, πœ‘π‘˜π‘₯ = 0, πœ‘π‘˜π‘¦ = 0, (π‘₯, 𝑦) ∈ πœ•π›Ί. 304 Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling Because of awkward form, governing equations of GCT are not presented. Detailed description of the theory can be found in [7]. It is necessary to note, that plate theories are approximations to elasticity theory and bring about their own errors which are important to estimate. 3 hp - Version of CLS Method Implementation aspects of CLS method are similar to the CM. Consider a general boundary value problem for a linear elliptic system in a rectangular domain 𝛺 = [0, π‘Ž] Γ— [0, 𝑏]: 𝐿𝑒(π‘₯, 𝑦) = 𝑓 (π‘₯, 𝑦), (π‘₯, 𝑦) ∈ 𝛺, 𝐿𝑏𝑛𝑑 𝑒(π‘₯, 𝑦) = 𝑔(π‘₯, 𝑦), (π‘₯, 𝑦) ∈ πœ•π›Ί. Let us define a grid with non-overlapping rectangular cells 𝛺 π‘˜ (π‘˜ = 1, . . . 𝐾). In each cell 𝛺 π‘˜ we introduce local variables (𝛼1π‘˜ , 𝛼2π‘˜ ), which are associated with global variables (π‘₯, 𝑦) in the Cartesian coordinate system by π‘₯ βˆ’ π‘₯*π‘˜ 𝑦 βˆ’ 𝑦 *π‘˜ 𝛼1π‘˜ = , 𝛼2π‘˜ = , π‘‘π‘˜1 π‘‘π‘˜2 where 2π‘‘π‘˜1 , 2π‘‘π‘˜2 β€” sizes of cell in π‘₯ and 𝑦 directions, (π‘₯*π‘˜ , 𝑦 *π‘˜ ) β€” the coordinates of the cell centers. Local variables are varying in canonical interval 𝛼1π‘˜ ∈ [βˆ’1, 1], 𝛼2π‘˜ ∈ [βˆ’1, 1]. The upper index π‘˜, that indicates the cell’s number, will be omitted further. In this version of CLS method approximate solution in the cell is represented in form of direct product of single variable basis functions : 1 βˆ’1 𝑁 π‘βˆ‘οΈ βˆ‘οΈ2 βˆ’1 𝑒(𝛼1 , 𝛼2 ) = 𝑐𝑖1 𝑖2 πœ‘π‘–1 (𝛼1 )πœ‘π‘–2 (𝛼2 ). (7) 𝑖1 =0 𝑖2 =0 Functions πœ‘π‘– are chosen as the Chebyshev polynomials of the first kind 𝑇𝑛 πœ‘π‘–1 (𝛼1 ) = 𝑇𝑖1 (𝛼1 ), πœ‘π‘–2 (𝛼2 ) = 𝑇𝑖2 (𝛼2 ). In previous paper [8] we used cardinal functions in Lagrange-like form 1 βˆ’1 π‘βˆοΈ 2 βˆ’1 π‘βˆοΈ 𝛼1 βˆ’ (𝛼1 )π‘π‘œπ‘™ π‘š 𝛼2 βˆ’ (𝛼2 )π‘π‘œπ‘™ 𝑙 πœ‘π‘–1 (𝛼1 ) = π‘π‘œπ‘™ βˆ’ (𝛼 )π‘π‘œπ‘™ , πœ‘π‘–2 (𝛼2 ) = π‘π‘œπ‘™ βˆ’ (𝛼 )π‘π‘œπ‘™ , (8) π‘š=0 (𝛼1 ) 𝑖1 1 π‘š 𝑙=0 (𝛼 2 )𝑖2 2 𝑙 π‘šΜΈ=𝑖1 𝑙̸=𝑖2 (οΈ€ )οΈ€ where (𝛼1 )π‘π‘œπ‘™ π‘π‘œπ‘™ π‘š , (𝛼2 )π‘š are local coordinates of collocation points. But in prac- tice polynomials (8) require a large number of arithmetic operations and lead to complex expressions when differentiating. In this sense, Chebyshev polynomials are more convenient choice. To determine the unknown coefficients in representation (7) for each cell let us write down the equation of three types 305 Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling (οΈ€ )οΈ€ – collocation equations at the collocation point 𝛼1π‘π‘œπ‘™ , 𝛼2π‘π‘œπ‘™ 𝐿𝑒(𝛼1π‘π‘œπ‘™ , 𝛼2π‘π‘œπ‘™ ) = 𝑓 (𝛼1π‘π‘œπ‘™ , 𝛼2π‘π‘œπ‘™ ); (9) (οΈ€ 𝑏𝑛𝑑 𝑏𝑛𝑑 )οΈ€ – boundary equations at given point 𝛼1 , 𝛼2 , on boundary πœ•π›Ί π‘˜ , adjacent to πœ•π›Ί 𝐿𝑏𝑛𝑑 𝑒(𝛼1𝑏𝑛𝑑 , 𝛼2𝑏𝑛𝑑 ) = 𝑔(𝛼1𝑏𝑛𝑑 , 𝛼2𝑏𝑛𝑑 ); (10) – matching conditions on interface between neighbour cells at given points (𝛼1π‘šπ‘Žπ‘‘ , 𝛼2π‘šπ‘Žπ‘‘ ) πΏπ‘šπ‘Žπ‘‘ 𝑒(𝛼1π‘šπ‘Žπ‘‘ , 𝛼2π‘šπ‘Žπ‘‘ ) = πΏπ‘šπ‘Žπ‘‘ π‘’π‘Žπ‘‘π‘— (𝛼1π‘šπ‘Žπ‘‘ , 𝛼2π‘šπ‘Žπ‘‘ ), (11) Here π‘’π‘Žπ‘‘π‘— – solution defined in neighbour cell 𝛺 π‘Žπ‘‘π‘— . Matching conditions πΏπ‘šπ‘Žπ‘‘ usually require the continuity of the solutions and the necessary number of its derivatives along the normal to the boundary of the cell. In this version of CLS method the local coordinates of collocation points are roots of Chebyshev polynomials (𝑖1 = 1, . . . , 𝑁1 , 𝑖2 = 1, . . . , 𝑁2 ) (οΈ€ )οΈ€ (𝛼1 )π‘π‘œπ‘™ π‘π‘œπ‘™ 𝑖1 , (𝛼2 )𝑖2 = (𝑑𝑖11 , 𝑑𝑖22 ), where 𝑑𝑖11 and 𝑑𝑖22 – roots of Chebyshev (οΈ€polynomials)οΈ€ of 𝑁1 and 𝑁2 degree re- spectively. By the same way we define 𝛼1𝑏𝑛𝑑 , 𝛼2𝑏𝑛𝑑 and (𝛼1π‘šπ‘Žπ‘‘ , 𝛼2π‘šπ‘Žπ‘‘ ) on cell boundaries (𝑖1 = 1, . . . , 𝑁1 , 𝑖2 = 1, . . . , 𝑁2 ) (βˆ’1, 𝑑𝑖22 ), (1, 𝑑𝑖22 ), (𝑑𝑖11 , βˆ’1), (𝑑𝑖11 , 1). Thus, for 𝑁1 𝑁2 unknown coefficients in cell we use 𝑁1 𝑁2 collocation equa- tions appended by equations on cell boundary. Thus, in hp - version of CLS method corresponding SLAE becomes overdetermined. In this version of method approximate solution form does not satisfy to boundary and matching conditions identically, so residual 𝑅(π‘₯) (2) must contain not only collocation equations, but boundary conditions (10) and matching conditions (11) too. In particular, this allows us to consider the boundary conditions in complex form. To solve the overdetermined linear systems in the least squares sense (2) we use 𝑄𝑅 factorization of its matrix, implemented by Householder method. In the case of large linear systems we use domain decomposition method [9]. This allows to reduce the solution of the problem in whole region to iterative process through subdomains with computational complexity is much smaller than the original problem for the region. For linear systems in a subdomain we use Householder method again. In this case special matching conditions between subdomains are used. For example, the continuity of function and its first derivative at the boundary of the cell 𝛺 π‘˜ can be written as πœ•π‘’ πœ•π‘’π‘Žπ‘‘π‘— 𝑒 + 𝑝1 = π‘’π‘Žπ‘‘π‘— + 𝑝1 , πœ•π‘› πœ•π‘› where 𝑒 is solution in the cell at the current iteration; π‘’π‘Žπ‘‘π‘— β€” solution in the neighbor cell; 𝑛 β€” the outer normal to the boundary 𝛺 π‘˜ . For plate theories, 306 Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling that may contain derivatives up to the 4th order, we can additionally require the continuity of a linear combination of the second and third derivatives of the approximate solution: πœ•2𝑒 πœ•3𝑒 πœ• 2 π‘’π‘Žπ‘‘π‘— πœ• 3 π‘’π‘Žπ‘‘π‘— 2 + 𝑝2 3 = 2 + 𝑝2 . πœ•π‘› πœ•π‘› πœ•π‘› πœ•π‘›3 The choice of weights 𝑝1 , 𝑝2 can affect the properties of the numerical solutions and speed of convergence of the iterative process. 4 Numerical Experiments In all numerical experiments a single cell that coincides with the entire domain is used. Further we will consider only square plates π‘Ž = 𝑏 = 1 m. To demonstrate the capabilities of hp - version of CLS method let us consider the problem with the known exact solution. The last equation of the system (6) is similar to the Kirchhoff-Love plate theory equation for bending of a homogeneous orthotropic plates πœ• 4 𝑀0 πœ• 4 𝑀0 πœ• 4 𝑀0 𝐷11 4 + (2𝐷12 + 4𝐷66 ) 2 2 + 𝐷22 = π‘ž0 . πœ•π‘₯ πœ•π‘₯ πœ•π‘¦ πœ•π‘¦ 4 Consider 3-ply simply-supported square laminate under uniform transverse load π‘ž0 with stiffness coefficients defined by (3). This problem can be solved by Fourier method [10]. In this case maximum deflection is observed in the center of the plate and if β„Ž = 0.01 m, then the deflection value for Fourier method solution (sum of first 2500 members) is 𝑀0 (0.5, 0.5) 8 𝑀* = 10 = 9.8577127. π‘ž0 Deflections in the center of the plate are calculated by hp – version of CLS method are shown in Table 1. Table 1. Deflections at the center of the plate obtained by the hp – version of CLS method, β„Ž = 0.01. 𝑁1 Γ— 𝑁2 𝑀0 (0.5, 0.5)/π‘ž0 Β· 108 10 Γ— 10 9.8579 20 Γ— 20 9.85773 30 Γ— 30 9.857715 Table 1 demonstrates high accuracy of numerical results, obtained by hp – version of CLS method even for the differential equation with derivatives of 307 Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling 4-th order of the unknown functions. Thus, hp – version of CLS method has no difficulties when working with differential equations containing derivatives of high orders like (3). Consider another formulation of the problem. Let us use the free edge con- ditions on the one of the edges (π‘₯ = 1) in the previous problem. πœ• 2 𝑀0 πœ• 2 𝑀0 𝐷11 + 𝐷 12 = 0, πœ•π‘₯2 πœ•π‘¦ 2 πœ• 3 𝑀0 πœ• 3 𝑀0 𝐷11 + (𝐷 12 + 2𝐷 66 ) = 0. πœ•π‘₯3 πœ•π‘₯2 πœ•π‘¦ These boundary conditions by the use of term (2) implements with no additional effort in the CLS method. Fig. 2 shows the deformed shape of the plate with the free edge. Fig. 2. Simply-supported plate with free edge. Now let us consider the bending of clamped 3-ply laminates with different relative thickness β„Ž/π‘Ž, as described in Section 2. Stress and displacement fields calculation for such plates will be carried out within framework of three theories described above. Fig. 3. Stress fields for absolute value of 𝜎π‘₯π‘₯ and πœŽπ‘¦π‘¦ in 3-ply laminate for β„Ž/π‘Ž=0.02. 308 Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling Brief analysis of solution shows that 𝜎π‘₯π‘₯ prevails in the stress state of outer layers (Fig. 3). And maximum absolute values are observed in the vicinity of the clamped edges (π‘₯ = 0, 1). Similar conclusion are true for the middle layer of the plate and component πœŽπ‘¦π‘¦ , which is associated with an orientation of transversely isotropic material. Further we will use the following normalized quantities π‘˜ π‘˜ π‘˜ 𝜎π‘₯π‘₯ πœŽπ‘¦π‘¦ 𝜎π‘₯π‘₯ 𝜎 π‘˜ Β―π‘₯π‘₯ = , 𝜎 π‘˜ ¯𝑦𝑦 = , 𝜎 π‘˜ Β―π‘₯𝑦 = 10βˆ’2 , π‘ž0 𝑆 2 π‘ž0 𝑆 2 π‘ž0 𝑆 2 π‘˜ 𝜎π‘₯π‘₯ 𝑧 𝑀 Β―0 = 10βˆ’9 , 𝑧¯ = . π‘ž0 β„Ž Table 2. Stresses and deflection in 3-ply laminate. The results of calculations carried out in the framework of the CLPT, FSDT and GCT plate theories. Sign (%) is used for relative percentage deviation from GCT. β„Ž/π‘Ž GCT FSDT CLPT FSDT (%) CLPT (%) 3 𝜎 Β―π‘₯π‘₯ (π‘Ž, 0, β„Ž/2) 0.1 0.975 0.475 0.579 51.2 40.6 0.05 0.569 0.545 0.579 4.28 1.61 0.02 0.576 0.574 0.578 0.30 0.42 0.01 0.578 0.578 0.578 0.05 0.13 2 𝜎 ¯𝑦𝑦 (0, π‘Ž, β„Ž/4) 0.1 0.636 0.662 0.456 4.21 28.2 0.05 0.605 0.536 0.456 11.5 24.6 0.02 0.483 0.471 0.456 2.65 5.72 0.01 0.463 0.460 0.455 0.64 1.82 3 𝜎 Β―π‘₯𝑦 (3/4π‘Ž, 3/4π‘Ž, β„Ž/2) 0.1 0.983 0.675 0.006 31.3 38.4 0.05 0.632 0.647 0.006 2.28 4.20 0.02 0.616 0.616 0.006 0.04 1.74 0.01 0.609 0.609 0.006 0.02 0.62 𝑀(π‘Ž/2, Β― π‘Ž/2, β„Ž/2) 0.1 –0.751 –0.603 –0.208 31.30 38.36 0.05 –2.650 –2.538 –1.667 2.28 4.20 0.02 –28.523 –28.335 –26.044 0.04 1.74 0.01 –213.312 –212.974 –208.354 0.02 0.62 Table 2 shows calculated deflections and stresses in vicinities of their maxi- mum absolute values. We assume GCT as the most accurate among considered theory, because its hypothesis is the most suitable in given structure [8, 7]. There- fore, we will treat it as a reference. 309 Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling Despite the simplicity, the CLPT may be used for very thin laminated plates. If the 5% deviations from GCL are admissible, CLPT can be applied for β„Ž/π‘Ž < 0.02 case. More accurate results are obtained within the framework of FSDT. It can be used for plates with aspect ratio β„Ž/π‘Ž < 0.05. For thicker plates GCT theory differs significantly. In presented calculations all functions are approximated with 𝑁1 = 𝑁2 = 17 that allowed to get three guaranteed digits for the maximum deflection. In this case, the CLPT is required to determine 17 Β· 17 Β· 3 = 867 unknown coefficients, FSDT – 17 Β· 17 Β· 5 = 1445, and GCT – 2601. It is interesting that for β„Ž/π‘Ž < 0.02 CLPT results do not differ significantly from those of FSDT. In particular, it means that for plates with β„Ž/π‘Ž < 0.02 CLPT is preferable, because of lower computational efforts. Fig. 4. Stresses distribution along 𝑧 coordinate 𝜎 Β―π‘₯π‘₯ (π‘Ž, 0, 𝑧¯), 𝜎 Β―π‘₯π‘₯ (π‘Ž/2, π‘Ž/2, 𝑧¯) and 𝜎 ¯𝑦𝑦 (0, 𝑏, 𝑧¯), 𝜎 ¯𝑦𝑦 (π‘Ž/2, π‘Ž/2, 𝑧¯) in 3-ply laminate for β„Ž/π‘Ž = 0.02. Fig. 4 shows the distribution of stresses 𝜎π‘₯π‘₯ and πœŽπ‘¦π‘¦ along the 𝑧 coordinate (thickness) for different points (π‘₯, 𝑦): in the center of the plate and on the border where they reach the maximum values. The absolute maximum stresses values are observed at the outer sutfaces of layers. And the absolute values at the edges of the plate exceed values at the center of the plate a few times, that is true for both the stress tensor components. Presented numerical results shows that hp - version of CLS method can be successfully applied to problems of mechanics of laminates anisotropic rectan- gular plates within framework the various plates theories. Term (2) allows us to 310 Mathematical and Information Technologies, MIT-2016 β€” Mathematical modeling consider a wide class of boundary value problems including complex boundary conditions. Moreover p - approach allows to obtain high accuracy of numerical solutions at low computational efforts. References 1. Russell, R.D., Shampine, L.F.: A collocation method for boundary value problems. Numer. Math. 19, 1–28 (1972) 2. Shapeev, V.P., Isaev, V.I., Idimeshev, S.V.: The collocations and least squares method:application to numerical solution of the Navier-Stokes equations. CD-ROM Proc. of the 6th ECCOMAS (Austria, Vienna, September 10–14, 2012). (2012) 3. Shapeev, V.P.: Collocation and Least Residuals Method and Its Applications. EPJ Web of Conferences 108, 01009, 1–12 (2016) 4. Bochev, P., Gunzburger, M.: Least Squares Finite Element Methods. Springer, Berlin (2009) 5. Pagano, N.J., Hatfield, H.J.: Elastic Behavior of Multilayered Bidirectional Com- posites. AIAA Journal. 10, 7, 931-933 (1972) 6. Reddy, J.N. Mechanics of Laminated Composite Plates and Shells. CRC Press (2003) 7. Grigoluk, E.I. Kulikov, T.M.: Multilayered Reinforced Shells. Design of Pneumatic Tires. Mashinostroyenie, Moscow (1988). 8. Golushko, S.K., Idimeshev, S.V., Shapeev, V.P.: Development and application of collocations and least residuals method to the solution of problems in mechanics of anisotropic laminated plates. Computational technologies. 19, 5, 24–36 (2014) 9. Golushko, S.K., Idimeshev, S.V., Shapeev, V.P.: Application of collocations and least residuals method to problems of the isotropic plates theory. Computational technologies. 18, 6. 31–43 (2013) 10. Timoshenko, S. and Woinowsky-Krieger, S.: Theory of plates and shells. McGraw– Hill, New York (1959) 311