On the formulation of a high-order discontinuous finite element method based on orthogonal polynomials for laminated plate structures

ARTICLE IN PRESS

JID: MS

[m5GeSdc;August 18, 2017;9:38]

International Journal of Mechanical Sciences 000 (2017) 1–19

Contents lists available at ScienceDirect

International Journal of Mechanical Sciences journal homepage: www.elsevier.com/locate/ijmecsci

On the formulation of a high-order discontinuous finite element method based on orthogonal polynomials for laminated plate structures Tianyu Li1,∗, Rakesh K. Kapania2 Kevin T. Crofton Department of Aerospace and Ocean Engineering, Virginia Polytechnic Institute and State University, Blacksburg, VA 24060, USA

a r t i c l e

i n f o

Keywords: Laminated plate Discontinuous finite element Diagonal mass matrix Mode shape function Consistent Orthogonal Basis Function Space Stress analysis

a b s t r a c t Laminated plate structures are analyzed by a discontinuous finite element method with emphasis on determining the transverse shear and normal stress components at the interface of adjacent layers accurately. A Consistent Orthogonal Basis Function Space is used for the interpolation of the displacement field and the traction field between two adjacent layers. The mass matrix of the laminated plate becomes diagonal. Moreover, it is observed that the basis functions are very similar to the vibration mode shapes, even through we do not solve any eigenvalue problem in their generation. These basis functions are uniquely determined by the structure’s configuration and associated boundary conditions. The stress field between the layers can be accurately calculated, even for the region near the boundaries that might have sharp stress gradients. Several numerical examples are studied with different boundary conditions. The results for both the deformation and the stress components are compared with the traditional finite element method, especially in terms of the number of degrees-of-freedom (DOF) used in the proposed method and the classic FEM. It is observed that the proposed method is able to use a much fewer number of DOF than that of commercial FEM software (ANSYS etc) to obtain accurate solutions to both the deformation of the plate and the stress field between adjacent layers. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction Laminated plates have wide applications due to their high strengthto-weight and stiffness-to-weight ratios. In the analysis of laminated structures, the finite element method has proven to be a highly reliable approach. The material on nonlinear finite element is available in [1]. The laminated beam/plate development history is described in [2,3]. Static and dynamic analysis with large deformation are both included. The constitutive theory used in the present paper is based on the material available in [4]. In [5], a laminated beam finite element is developed wherein the shear force effect is considered. The laminated plate buckling problem is solved in [6], wherein the geometrical imperfections are considered. The layerwise theory is applied to solve the laminated plate vibration problem in [7]. The stiffened panel buckling analysis is available in [8]. The global-local discontinuous shell element for large strains and rotation is developed in [9], where the shear deformation is also considered. In [10,11], a multi-variable Hellinger-Reissner variation principle is applied for laminated plate structures and the shear force is used as a degree-of-freedom. A layer-wise theory is used for studying displacements in the laminated plates in [12]. Laminated plate

∗

1 2

analysis with continuous shear force across the thickness direction is available in [13]. In [14], the buckling and vibration of a laminated plate by using different plate theories are presented. The discontinuous layerwise displacement model is applied to analyze the composite shells in [15]. In [16], the Reissner-mixed variational theorem is applied to analyze the static and free vibration response of cross-ply laminated plates, where the stresses and displacement are unknowns. In [17], an enriched shell element based on floating node method is applied to analyze the delamination of laminated plate, where discontinuous displacement is also considered. In [18], the impact problem of laminated plate is solved, with emphasis on the study of sensitivity of material parameters. The localization of buckling modes to laminated plates is studied in [19] by using Mindlin plate finite element. In [20], the kp-Ritz method is applied to analyze the laminated CNT reinforced plates, where the high-order displacement is also employed. A high-order hybrid discontinuous finite element method is developed here for the laminated plate analysis with large deformation and finite strain. The proposed method is very accurate for the laminated plates while using a much fewer number of DOF (degree-of-freedom) when comparing with the traditional FEM. The discontinuous finite element is used for the laminated plate because:

Corresponding author. E-mail address: [email protected] (T. Li). Graduate student Norris and Wendy Mitchell Endowed Professor

http://dx.doi.org/10.1016/j.ijmecsci.2017.08.006 Received 19 March 2017; Received in revised form 31 July 2017; Accepted 1 August 2017 Available online xxx 0020-7403/© 2017 Elsevier Ltd. All rights reserved.

Please cite this article as: T. Li, R.K. Kapania, International Journal of Mechanical Sciences (2017), http://dx.doi.org/10.1016/j.ijmecsci.2017.08. 006

JID: MS

ARTICLE IN PRESS

T. Li, R.K. Kapania

[m5GeSdc;August 18, 2017;9:38] International Journal of Mechanical Sciences 000 (2017) 1–19

(1) The stress field between adjacent layers is accurately solved; (2) The discontinuous finite element allows discontinuous basis functions, which provides the freedom to further optimize the basis function space; For the discontinuous finite element in this paper, the basis function space of displacement is further developed based on orthogonal polynomial theory. An accurate diagonal mass matrix is obtained without any numerical approximations (such as lumped mass matrix or the artificial diagonal mass matrix by Gauss-Lobatto numerical integration). The lumped diagonal mass matrix is studied in [21]. The spectral element method with Legendre polynomials is studied in [22] with a diagonal mass matrix. The dynamical study of the Legendre spectral element is presented in [23]. In reference [24], the Legendre spectral element method is applied to analyze composite plates. In [25], incompressible continua is studied with using Legendre spectral element. The spectral element method for composite strips is studied in [26]. The diagonal mass matrix is obtained for all the references [21–26] for solving various problems in solid mechanics. However, the uniqueness of the structure configuration and the boundary conditions are not considered in the development of the basis functions in previous work. They use the same basis functions for all structures. In this paper, for a given structure and the associated boundary conditions, a unique basis function space is generated. For different structures or boundary conditions, the basis function spaces will be different. Additional references of analytical and numerical methods as applied to the analysis of laminated plates are presented here. Analytical solutions for laminated plates are available in [27,28], in which the analytical solutions for static and dynamical analysis of sandwich plates are developed. Higher order refined theories are applied to model the sandwich behavior. Shear deformable plate theories are widely used in analyzing laminated plates. In [29], the finite element method based on higher order shear deformable theory is applied to analyze sandwich plates. In [30], a simple higher-order theory is used for analyzing laminated composite plates, in which the linear shear deformable plate theory and the parabolic shear strain distribution are both included. In [31], a third order shear deformable theory for functionally graded plates is developed. In [32], a postbuckling analysis of carbon nanotube reinforced functionally graded plates is presented. The so-called Zig-zag theories [33] are also widely used in laminated plates. In a Zig-zag theory, the displacement basis functions are designed such that the compatibility condition of shear strains are satisfied at the interfaces. The piecewise linear functions are used in the thickness direction. For a through study of the mechanics on laminated plates, a reader is referred to [34]. For more papers on laminated plates, it is referred to the review papers [35,36]. In this paper, a method is proposed based on 3-D continua mechanics theory. In the thickness direction, the discontinuous finite element method is applied to satisfy the compatibility condition. The interface stress field is also solved. The basis functions based on orthogonal polynomials (Consistent Orthogonal Basis Functions) are applied to analyze laminated plates. The Consistent Orthogonal Basis Function Space has the following advantages:

Fig. 1. The notations of laminated plate structure.

mathematical functions, but they also have a physical meaning. In this paper, the Consistent Orthogonal Basis Function Space is discussed especially for the laminated plate structure. We will deal with several cases of boundary conditions in Section 6. Different from those mixed finite element methods that use the displacement and the stress as DOFs, for the proposed discontinuous finite element method, only the stresses on the interface are solved. This choice is suitable for the laminated plates, in which case the interface stresses are more important than the stresses elsewhere. The basis functions of the interface stresses are also based on the Consistent Orthogonal Basis Functions developed in this paper. The proposed method has the potential to be generalized to solve the delamination of laminated plates. This will be done in future. The outline of the paper is as follows. In Section 2, the kinematics of the laminated plate structure is presented. In Section 3, the weak form of the discontinuous finite element method is discussed. In Section 4, the Consistent Orthogonal Basis Function Space is developed. The diagonal mass matrix is obtained and the generation procedure of the basis functions is given. Four specific examples are also presented. In Section 5, the nonlinear finite element implementation details are given. The second Piola-Kirchhoff stress and the Green strain are used in the calculation of strain energy. The constitutive relation is set up between Jaumann stress and linear strain rate. In Section 6, the numerical tests are presented. In Section 7, the conclusion from this study is given. 2. Kinematics Assume the laminated plate structure has nl layers. The top surface and bottom surface are denoted as Γ± , 𝑘 = 1, 2, 3, ⋯ , 𝑛𝑙 . We also define 𝑡,𝑘 Γ𝑖𝑡,𝑘 as the interface between the kth layer and the k + 1th layer. Thus, we have Γ𝑖𝑡,𝑘 = Γ+ = Γ− . The space at time t + dt for the kth layer is 𝑡,𝑘 𝑡,𝑘+1 denoted as Ωt,k , as Fig. 1 shows. For the laminated plate structure, we assume that the natural coordinate in the thickness direction is defined locally for every layer:

(1) It generates the basis functions that are quite similar to the mode shapes of the structure under consideration, yet these are obtained without solving the eigenvalue problem for the structure; (2) The use of the basis function space leads to a very high accuracy by using only a few number of DOFs; (3) A diagonal mass matrix is accurately obtained, which makes explicit dynamic analysis faster because of total decoupling of the degree-of-freedom.

𝜉3𝑘 =

𝑧 − 𝑧𝑘 ∈ [0, 1] 𝑧𝑘+1 − 𝑧𝑘

(2-1)

For the plate dimension, we have: 𝑥 𝑎 𝑦 𝜉2 = 𝑏 𝜉1 =

(2-2a,b)

where a, b are the length and width of the laminated plate, and x, y, z are the global coordinates, zk ,k = 0, 1, ⋅⋅⋅, nl is the z coordinate for every layer.

With the development of basis functions similar to the mode shapes, a faster convergent rate is obtained. The basis functions are not only 2

ARTICLE IN PRESS

JID: MS

[m5GeSdc;August 18, 2017;9:38]

T. Li, R.K. Kapania

International Journal of Mechanical Sciences 000 (2017) 1–19

The displacement trial function spaces are next discussed. The function space is defined as below; for i = 1, 2, 3: 𝑚

𝑚 𝑉 = {𝜙𝑢𝑖0 |𝜙𝑚 𝑢 (𝜉) ∈ 𝑃 (𝜉), 𝜉 ∈ [0, 1], 𝑚 = 1, 2, 3, ⋯ , 𝑚0 } 𝑖

traction field is also a primary variable, which makes the present approach to be different from the traditional displacement finite element method. The advantage of the current approach is that the traction field is continuous at the interfaces, even for the large deformation analysis. For the traditional finite element method, the stress obtained from the constitutive relations is generally not continuous at the interfaces. For the kth layer, the discontinuous finite element weak form will be:

(2–3)

where Pl is the polynomial function space whose highest order term is l. Then, a tensor product is used to give the multi-dimensional displacement field as: 𝜉

𝜉

𝜉

𝑢𝑘1 = 𝑅0 𝑓𝑢11 𝑓𝑢12 𝑓𝑢13 𝜉

𝜉

𝜉

𝑢𝑘2 = 𝑅0 𝑓𝑢21 𝑓𝑢22 𝑓𝑢23 𝜉

𝑚0 𝑛0 𝑠0 ∑ ∑∑ 𝑚=1 𝑛=1 𝑠=1 𝑚0 𝑛0 𝑠0

∑ ∑∑

𝑚=1 𝑛=1 𝑠=1 𝑚0 𝑛0 𝑠0

∑ ∑∑

𝜉

𝑢𝑘3 = 𝑅0 𝑓𝑢31 𝑓𝑢32 𝑓𝑢𝜉

3

𝑚=1 𝑛=1 𝑠=1

𝑛𝑙 ∑

( ) 𝑛 ( ) 𝑠 ( 𝑘) 𝐴𝑚𝑛𝑠 𝜙𝑚 𝑢 𝜉1 𝜓𝑢 𝜉2 𝜑𝑢 𝜉3 1

1

𝑘=1

1

−

( ) 𝑛 ( ) 𝑠 ( 𝑘) 𝐵𝑚𝑛𝑠 𝜙𝑚 𝑢 𝜉1 𝜓𝑢 𝜉2 𝜑𝑢 𝜉3 2

2

𝑛𝑙 ⟨ ∑ 𝑘=1

⟩ 𝑘

𝛿𝑢𝑘𝑖 , 𝑇̄𝑖

0,𝑘

Γ𝜎𝑡,𝑘

−

𝑘=1 𝑛𝑙 ⟨

∑

𝑘=1

⟩ 𝑘

𝛿𝑢𝑘𝑖 , 𝑇𝑖

Ω𝑡,𝑘

Γ± 𝑡,𝑘

−

𝑛𝑙 ⟨ ∑ 𝑘=1

−

𝑛𝑙 ∑ 𝑘=1

3

3

⟨ ⟩ 𝑘 𝛿 𝑢̃ 𝑘𝑇 𝑤𝑘𝑇 𝑖 , 𝑏𝑖 Ω

𝑡,𝑘

⟩ 𝑘

𝛿𝑇𝑖𝑘 , 𝑢𝑖

2

( ) 𝑛 ( ) 𝑠 ( 𝑘) 𝐶𝑚𝑛𝑠 𝜙𝑚 𝑢 𝜉1 𝜓𝑢 𝜉2 𝜑𝑢 𝜉3 3

⟩ 𝑛𝑙 ⟨ ⟨ ⟩ ∑ 𝑘 𝑘 𝛿 𝑢̃ 𝑘𝑇 𝑤𝑘𝑇 𝛿𝑡𝑡+𝑑𝑡 𝜀𝑘𝑖𝑗 ,𝑡𝑡+𝑑𝑡 𝑆𝑖𝑗𝑘 𝑖 , 𝜌0 𝑢̈ 𝑖 Ω +

Γ± 𝑡,𝑘

=0 (3-1)

where 𝜌0 is the density field in the initial configuration space, 𝑤𝑘𝑖 is the weighting function for the kth layer in the ith dimension, 𝑡𝑡+𝑑𝑡 𝜀𝑘𝑖𝑗

(2–4a,b,c)

is the Green strain field for the kth layer, 𝑡𝑡+𝑑𝑡 𝑆𝑖𝑗𝑘 is the corresponding second Piola-Kirchhoff stress tensor, 𝑏𝑘𝑖 is the body force component in the ith dimension for the kth layer, 𝑇̄𝑖𝑘 is the prescribed traction field component on the stress boundary Γ𝜎𝑡,𝑘 for the ith layer, 𝑇𝑖𝑘 is the traction field defined on the interface to be determined as a primary variable, 𝛿 𝑢̃ 𝑖 is the variant of the displacement DOF vector 𝑢̃ 𝑖 .

where 𝑢𝑘1 , 𝑢𝑘2 , 𝑢𝑘3 are denoted as the x, y, z displacement; 𝜉1 , 𝜉2 , 𝜉3𝑘 are the natural coordinates; m0 ,n0 ,s0 are the highest orders for 𝜉1 , 𝜉2 , 𝜉3𝑘 ,respectively; R0 is a reference length with dimension of length, which is only a constant used to non-dimensionalize the degrees of 𝑘 , 𝐶 𝑘 ; 𝜙𝑚 (𝜉 ) ∈ 𝑉 , 𝑖 = 1, 2, 3 are the polynomial funcfreedom 𝐴𝑘𝑚𝑛𝑠 , 𝐵𝑚𝑛𝑠 𝑚𝑛𝑠 𝑢𝑖 1 tions about the natural coordinate 𝜉 1 for ui ; 𝜓𝑢𝑛 (𝜉1 ) ∈ 𝑉 , 𝑖 = 1, 2, 3 and 𝑖

𝜑𝑠𝑢 (𝜉3𝑘 ) ∈ 𝑉 , 𝑖 = 1, 2, 3 are the displacement basis functions about the nat-

4. Consistent Orthogonal Basis Function Space

𝑖

ural coordinate 𝜉 2 and natural coordinate 𝜉3𝑘 , respectively. The specific 𝜉

formulation of 𝜙𝑢𝑖 , 𝜓𝑢𝑖 , 𝜑𝑢𝑖 will be discussed later in Section 4; 𝑓𝑢𝑖𝑗 is a displacement boundary condition function, which is used to define the displacement boundary condition about the natural coordinate 𝜉 j . If there is no displacement boundary condition about the natural coordi-

In this Section, the development of the Consistent Orthogonal Basis Function Space is presented [37]. For the sake of completeness of the mathematical formulation, the basic property of orthogonal polynomial is briefly reviewed. The orthogonal polynomials have the property that:

𝜉

nate 𝜉 j for ui ,𝑓𝑢𝑖𝑗 is defined as 1. For the laminated plate structure, the determinant of the Jacobian matrix is: det (𝐽 ) = 𝑡𝑘 𝑎𝑏 𝑡𝑘

⟨

3

𝑝0 = 1 𝑝 1 = 𝜉 − 𝛼0 ( ) 𝑝𝑛+1 = 𝜉 − 𝛼𝑛 𝑝𝑛 − 𝛽𝑛 𝑝𝑛−1 , 𝑛 ≥ 1

The interpolation of interface traction is given as: 𝜉

𝑇1𝑘 = 𝑇0 𝑓𝑇 1 𝑓𝑇 2 𝜉

1

𝜉

𝑇2𝑘 = 𝑇0 𝑓𝑇 1 𝑓𝑇 2 2

𝑇3𝑘

=

2

𝜉 𝜉 𝑇0 𝑓𝑇 1 𝑓𝑇 2 3 3

𝑚=1 𝑛=1 𝑚0 𝑛0

∑∑

𝑚=1 𝑛=1 𝑚0 𝑛0

∑∑

𝑚=1 𝑛=1

( ) 𝑛( ) 𝐴̃ 𝑚𝑛 𝜙𝑚 𝑇 𝜉1 𝜓𝑇 𝜉2 1

𝛼𝑛 =

( ) 𝑛( ) 𝐶̃𝑚𝑛 𝜙𝑚 𝑇 𝜉1 𝜓𝑇 𝜉2

𝛽𝑛 =

3

(2–7a,b,c)

𝑖

‖𝑝 ‖2 ‖ 𝑛 ‖𝑤

,𝑛 ≥ 0

‖𝑝 ‖2 ‖ 𝑛 ‖𝑤 ,𝑛 ≥ 1 ‖𝑝 ‖2 ‖ 𝑛−1 ‖𝑤

(4-3a,b)

In this case, an orthogonal polynomial series about coordinate 𝜉 is calculated. Once the weighting function is given, we can uniquely determine an orthogonal polynomial series in (4–1, 4–2). The displacement vector defined in (2–4a, b, c) is re-written in matrix form as:

where T0 is a reference traction used to non-dimensionalize the DOF 𝐴̃ 𝑚𝑛 , 𝐵̃ 𝑚𝑛 , 𝐶̃𝑚𝑛 , 𝜙𝑚 , 𝜓𝑇𝑛 ∈ 𝑉 are the basis functions of 𝜉 1 , 𝜉 2 for 𝑇𝑖 , 𝑖 = 𝑇 𝜉

⟨𝑥𝑝𝑛 , 𝑝𝑛 ⟩𝑤

2

3

(4-2a,b,c)

where

1

( ) 𝑛( ) 𝐵̃ 𝑚𝑛 𝜙𝑚 𝑇 𝜉1 𝜓𝑇 𝜉2 2

(4-1a,b,c)

where pi , pj are the i, j term of the orthogonal polynomial, respectively; w is the weighting function. It is noted that the weighting function w refers to the orthogonal polynomial only. On the other hand, once a weighting function w is given, we can always find an orthogonal polynomial series, such that (4–1a) holds.

3

𝑚0 𝑛0 ∑ ∑

= 𝛿𝑖𝑗

𝛿𝑖𝑗 = 0, 𝑖 ≠ 𝑗

𝑧𝑘

3

𝑤

𝛿𝑖𝑗 = 1, 𝑖 = 𝑗

Where = − is the thickness of the kth layer and the Jacobian matrix is defined as: 𝜕𝑦 𝜕𝑧 ⎤ ⎡ 𝜕𝑥 𝜕𝜉1 𝜕𝜉1 ⎥ ⎢ 𝜕𝜉1 ⎢ ⎥ 𝜕𝑦 𝜕𝑥 𝜕𝑧 ⎥ 𝐽 =⎢ (2-6) 𝜕𝜉2 𝜕𝜉2 ⎥ ⎢ 𝜕𝜉2 ⎢ ⎥ 𝜕𝑦 𝜕𝑧 ⎥ ⎢ 𝜕𝑥 ⎣ 𝜕𝜉 𝜕𝜉 𝜕𝜉 ⎦

1

⟩

(2-5)

𝑧𝑘+1

𝜉

𝑝𝑖 , 𝑝𝑗

𝑖

1, 2, 3, 𝑓𝑇 𝑗 is the traction boundary condition function of 𝜉𝑗 , 𝑗 = 1, 2 for 𝑖 𝑇𝑖 , 𝑖 = 1, 2, 3. It is noted that the traction expression in (2–7a, b, c) is defined on the interface plane 𝜉 1 -𝜉 2 .

𝑢 = 𝑁𝑢 𝑢̃

(4-4a)

where

3. Weak form

[ ] 𝑁𝑢 = 𝑁111 𝑁211 ⋯ 𝑁𝑚11 𝑁121 ⋯ 𝑁𝑚21 ⋯ 𝑁𝑚𝑛1 ⋯ 𝑁𝑚𝑛𝑘 [ ] 𝑢̃ = 𝑢̃ 𝑇111 𝑢̃ 𝑇211 ⋯ 𝑢̃ 𝑇𝑚11 𝑢̃ 𝑇121 𝑢̃ 𝑇121 ⋯ 𝑢̃ 𝑇𝑚21 ⋯ 𝑢̃ 𝑇𝑚𝑛1 ⋯ 𝑢̃ 𝑇𝑚𝑛𝑘

For the present method, the traction field Ti defined on the interface between adjacent layers are also introduced as a set of unknown variables to be interpolated and thus determined independently. The

(4-4b,c) 3

ARTICLE IN PRESS

JID: MS

[m5GeSdc;August 18, 2017;9:38]

T. Li, R.K. Kapania

International Journal of Mechanical Sciences 000 (2017) 1–19

In (4–4b), the matrix Nmnk is defined as:

where

𝑁𝑚𝑛𝑘 = 𝑁𝑚𝜙 𝑁𝑛𝜓 𝑁𝑠𝜑

(4-5a)

𝑚𝑎 =

where 𝑁𝑚𝜙

⎡𝑓𝑢𝜉1 𝜙𝑚 𝑢 ⎢ 1 1 =⎢ 0 ⎢ 0 ⎣

𝑁𝑠𝜑

𝑢̃ 𝑚𝑛𝑘

0

𝜉

𝑓𝑢21 𝜙𝑚 𝑢

2

0

⎡𝑓𝑢𝜉2 𝜓𝑢𝑛 ⎢ 1 1 =⎢ 0 ⎢ 0 ⎣

𝑁𝑛𝜓

⎤ ⎥ 0 ⎥ 𝜉1 𝑚 ⎥ 𝑓𝑢3 𝜙𝑢 ⎦ 3

0

𝜉

2

0

⎡𝑓𝑢𝜉3 𝜑𝑠𝑢 1 ⎢ =⎢ 0 ⎢ 0 ⎣

𝑚𝑏 =

𝑚𝑐 =

0

𝜉

𝑓𝑢23 𝜑𝑠𝑢

2

0

(4-5b,c,d)

[(

1

[(

1

[(

∫0

𝜉

𝑓𝑢21 𝜉

𝑓𝑢23 𝜉

𝑓𝑢31 𝜉

𝑓𝑢33

𝜉

∫

𝑁𝑢𝑇 𝑁𝑢 det (𝐽0 )𝑑𝜉1 𝑑𝜉2 𝑑𝜉3

∫0

(4–7)

𝑀1𝑛𝑑 ⎤ 𝑀2𝑛𝑑 ⎥⎥ ⋮ ⎥ 𝑀𝑛𝑑 𝑛𝑑 ⎥⎦

1

2

1

𝜉

det (𝐽0 ) = 𝑓det1 (𝜉1 )𝑓det2 (𝜉2 )𝑓det3 (𝜉3 )

]

𝜉

𝑘

∫0

𝜉

𝑓𝑢32

] 𝜉 𝑛 𝑛 𝑓det2 𝜓𝑢31 𝜓𝑢32 𝑑 𝜉2

𝑘

𝑓det3 𝜙𝑢31 𝜙𝑢32 𝑑 𝜉3

(4-15a,b,c)

𝜉

𝜉

𝜉

𝜉

)2

𝑓𝑢𝑖3

] 𝜉 𝑛 𝑛 𝑓det2 𝜓𝑢𝑖1 𝜓𝑢𝑖2 𝑑 𝜉2 = 𝛿𝑛1 𝑛2 , (𝑖 = 1, 2, 3)

)2

] 𝜉 𝑘 𝑘 𝑓det3 𝜙𝑢𝑖1 𝜙𝑢𝑖2 𝑑 𝜉3 = 𝛿𝑘1 𝑘2 , (𝑖 = 1, 2, 3)

(4–17)

𝜉

𝜉

𝜉

as the weight functions. Thus, for different structures (differ𝜉

𝑚0 ∑ 𝑚=1

𝐴𝑚 𝜙𝑚 𝑤 (𝜉1 )

𝜉

1

(4-16a,b,c)

𝜉

(𝑁𝑛𝜓 𝑁𝑛𝜓 𝑓det2 )𝑑𝜉2 1

2

1

∫0

𝜉

(𝑁𝑘𝜑 𝑁𝑘𝜑 𝑓det3 )𝑑𝜉3 1

(4–18)

will be defined as 𝑓𝑤1 = 𝜉12 (1 − 𝜉1 )2 , which means that the deflection is zero for the beam’s end points (𝜉1 = 0, 1). Based on the orthogonal basis function space developed in (4–15, 4–16), we have: 𝜙𝑚 𝑤 (𝜉1 ) to be orthogonal polynomials with the weight function 2 (𝜉1 (1 − 𝜉1 )2 )2 Then, we plot the first 6 basis functions of 𝜙𝑚 𝑤 (𝜉1 ) in Fig. 2. We see that the 𝜙𝑚 𝑤 (𝜉1 ) are almost identical with the mode shape functions of a fixed-fixed beam. If we change the boundary conditions to the fixed-free case, based on (4–15, 4–16), it is observed that 𝜙𝑚 𝑤 (𝜉1 ) are the orthogonal polynomials

2

(4–13) Substitute the (4–5a, b, c, d) into (4–13), thus, 0⎤ 0⎥ ⎥ 𝑚𝑐 ⎦

] 𝜉 𝑚 𝑚 𝑓det1 𝜑𝑢31 𝜑𝑢32 𝑑 𝜉1

where 𝜉1 = 𝐿𝑥 , L is the length of the beam. Am is the DOF to be solved, 𝜙𝑚 𝑤 (𝜉1 ) is the basis functions to be determined, the boundary conditions

(4–12)

𝜉

0 𝑚𝑏 0

[(

𝜉

where 𝑓det1 = 𝑎, 𝑓det2 = 𝑏, 𝑓det3 = 𝑡𝑘 Then, we have

⎡𝑚 𝑎 𝑀𝑖1 𝑖2 = 𝜌𝑠 𝑅20 ⎢ 0 ⎢ ⎣0

1

𝜉

𝑓𝑢𝑖2

𝑤 = 𝑓𝑤1

2

as:

∫0

)2

𝑘

𝜉

From (2–5 ), we can rewrite the determinant of the Jacobian matrix

2

𝑘

𝑓det3 𝜙𝑢21 𝜙𝑢22 𝑑 𝜉3

𝑓𝑢𝑖𝑗 , (𝑖 = 1, 2, 3, 𝑗 = 1, 2, 3)), there exists a unique basis function space, which is totally different from other structures. To illustrate the concepts, we give some simple examples, for which the closed-form representation of the mode shapes is available. We develop the orthogonal basis function space for each of these simple examples. It is very interesting that the orthogonal basis functions developed in this paper are quite identical to the closed-form mode shapes. For a one-dimensional uniform beam with fixed-fixed boundary conditions, the mid-surface displacement of deflection is assumed to be:

(4-10a,b)

1

1

𝜉

∫0

ent in the Jacobian matrix’s determinant 𝑓det𝑖 , 𝑖 = 1, 2, 3) and/or different boundary conditions (different in the boundary condition functions

where m1 , m2 , n1 , n2 , k1 , k2 are sub-script, and m0 , n0 , k0 are the order of basis functions, 𝑁𝑚𝜙1 𝑁𝑚𝜙2 𝑁𝑛𝜓1 𝑁𝑛𝜓2 𝑁𝑘𝜑 𝑁𝑘𝜑 are defined in (4–5b, c, d).

𝜉

[(

𝜉 𝜉 𝑓det3 (𝑓𝑢𝑖3 )2

(4-11a,b,c)

(𝑁𝑚𝜙 𝑁𝑚𝜙 𝑓det1 )𝑑𝜉1

[(

] 𝑛 𝜉 𝑛 𝑓det2 𝜓𝑢21 𝜓𝑢22 𝑑 𝜉2

𝜉

2

1 < 𝑘1 , 𝑘2 < 𝑘0

𝜉

]

𝜉

𝑓𝑢22

is easy to obtain in (4–2, 4–3) by taking 𝑓det1 (𝑓𝑢𝑖1 )2 , 𝑓det2 (𝑓𝑢𝑖2 )2 and

(4–9)

1 < 𝑛1 , 𝑛2 < 𝑛0

𝜉

] 𝜉 𝑚 𝑚 𝑓det1 𝜑𝑢21 𝜑𝑢22 𝑑 𝜉1

𝑀𝑖1 𝑖2 = 𝜌𝑠 𝛿𝑚1 𝑚2 𝛿𝑛1 𝑛2 𝛿𝑘1 𝑘2 𝐼

1 < 𝑚1 , 𝑚2 < 𝑚0

𝜉

1

)2

𝑘

In this case, the mass matrix is a diagonal matrix. For the specific formulations of 𝜙𝑢𝑖 (𝑖 = 1, 2, 3), 𝜓𝑢𝑖 (𝑖 = 1, 2, 3) and 𝜑𝑢𝑖 (𝑖 = 1, 2, 3), it

( ) ( ) 𝑖1 = 𝑚1 + 𝑛1 − 1 𝑚0 + 𝑘1 − 1 𝑚0 𝑛0 ( ) 0 ( ) 𝑖2 = 𝑚2 + 𝑛2 − 1 𝑚 + 𝑘2 − 1 𝑚0 𝑛0

𝜉

[(

] 𝜉 𝑛 𝑛 𝑓det2 𝜓𝑢11 𝜓𝑢12 𝑑 𝜉2

Thus (4–8)

(𝑁𝑚𝜙 𝑁𝑚𝜙 𝑁𝑛𝜓 𝑁𝑛𝜓 𝑁𝑘𝜑 𝑁𝑘𝜑 ) det (𝐽0 )𝑑𝜉1 𝑑𝜉2 𝑑𝜉3 2

1

∫0

where 1

)2

1

)2

[( ) ] 𝜉 2 𝜉 𝑚 𝑚 𝑓𝑢𝑖1 𝑓det1 𝜑𝑢𝑖1 𝜑𝑢𝑖2 𝑑 𝜉1 = 𝛿𝑚1 𝑚2 , (𝑖 = 1, 2, 3)

∫0

where J0 is the Jacobian matrix calculated from configuration and det (𝐽0 ) is the determinant of the Jacobian matrix J0 . Substitute (4–5a) into (4–7), the mass matrix can be written as:

∫ ∫ ∫

)2

𝑘

∫0

𝜉

𝑓𝑢12

𝑓det3 𝜙𝑢11 𝜙𝑢12 𝑑 𝜉3

𝜉

= 𝜌𝑠

⋯ ⋯ ⋱ ⋯

)2

𝜉

[(

𝜉

1

𝑀21 𝑀22 ⋮ 𝑀𝑛𝑑 2

)2

]

1

weighting function 𝑓det3 (𝑓𝑢𝑖3 )2 in the interval [0, 1], then

𝑁𝑢𝑇 𝑁𝑢 𝑑Ω

⎡ 𝑀11 ⎢𝑀 𝑀𝑚 = ⎢ 21 ⎢ ⋮ ⎢𝑀 ⎣ 𝑛𝑑 1

)2

] 𝜉 𝑚 𝑚 𝑓det1 𝜑𝑢11 𝜑𝑢12 𝑑 𝜉1

as the orthogonal series with weighting function 𝑓det2 (𝑓𝑢𝑖2 )2 in the interval [0, 1] and 𝜑𝑢𝑖 (𝑖 = 1, 2, 3) is selected as the orthogonal series with

∫

∫0

1

∫0

𝜉

𝑓𝑢13

)2

function 𝑓det1 (𝑓𝑢𝑖1 )2 in the interval [0, 1] and 𝜓𝑢𝑖 (𝑖 = 1, 2, 3) is selected

(4–6)

𝑀 𝑚 = 𝜌𝑠

1

[(

∫0

The mass matrix is defined as:

𝑀 𝑖 1 𝑖 2 = 𝜌𝑠

1

𝜉

𝑓𝑢11

If 𝜙𝑢𝑖 (𝑖 = 1, 2, 3) is selected as the orthogonal series with weighting

⎡𝐴𝑚𝑛𝑘 ⎤ = ⎢𝐵𝑚𝑛𝑘 ⎥ ⎢ ⎥ ⎣𝐶𝑚𝑛𝑘 ⎦

𝑀 𝑖 1 𝑖 2 = 𝜌𝑠

[(

∫0

⎤ ⎥ 0 ⎥ 𝜉3 𝑠 ⎥ 𝑓𝑢3 𝜑𝑢 ⎦ 3

0

1

∫0

0

𝑓𝑢22 𝜓𝑢𝑛

[(

∫0

⎤ ⎥ 0 ⎥ 𝜉 𝑓𝑢32 𝜓𝑢𝑛 ⎥⎦ 3

0

1

(4–14)

4

JID: MS

ARTICLE IN PRESS

T. Li, R.K. Kapania

[m5GeSdc;August 18, 2017;9:38] International Journal of Mechanical Sciences 000 (2017) 1–19

Fig. 2. Orthogonal basis functions for fixed-fixed uniform beam.

Fig. 3. Orthogonal basis functions for fixed-free uniform beam.

with the weight function (𝜉12 )2 . Thus, the plots of polynomials 𝜙𝑚 𝑤 (𝜉1 ) are presented in Fig. 3. From Fig. 3, we can tell that 𝜙𝑚 𝑤 (𝜉1 ) are quite identical to the mode shapes of a fixed-free beam. Let us analyze another example. In this case, the beam’s cross-section area is a function of axial coordinate: 𝐴(𝑥) = 𝐴0 (1 + 𝐿𝑥 ) Where we just assume that 𝐴0 = 1, 𝐿 = 1. In this case, the beam’s stiffness will be a function of x. Stiffness will increase from left to right. The fixed-fixed boundary conditions are applied. Thus, based on the 2 2 2 (4–15, 4–16), the weight of 𝜙𝑚 𝑤 (𝜉1 ) will be (𝜉1 (1 − 𝜉1 ) ) (1 + 𝜉1 ). Thus, the plots of polynomials 𝜙𝑚 ( 𝜉 ) are presented in Fig. 4 . 𝑤 1 From Fig. 4, it is interesting to observe that, these orthogonal basis functions have a larger value on the left part. For this beam, the left part’s stiffness is smaller than that of the right part (recall 𝐴(𝑥) = 𝐴0 (1 + 𝐿𝑥 ). This is valid with the mode shape functions, which will have a relatively larger magnitude for the relatively weaker part of the structure. If we change the boundary conditions to fixed-free, based on (4–15, 4–16), it is obtained that the functions 𝜙𝑚 𝑤 (𝜉1 )’s weight function is (𝜉12 )2 (1 + 𝜉1 ). In this case, the plots of polynomials 𝜙𝑚 𝑤 (𝜉1 ) are presented in Fig. 5. Comparing the plots in Fig. 5 with those in Fig. 3, it is found that the basis functions’ value is reduced on the right side. This is quite consistent

with the mode shapes because the beam’s stiffness on the right side is relatively higher. Thus, it is valid to state that, the orthogonal basis function space developed in this paper is a good representation of the mode shape functions. This is because that the basis functions are uniquely determined by the boundary conditions and the Jacobian matrix’s determinant. With the basis functions, the mass matrix is diagonal without any approximations. It is well known that the mode shape functions are also dependent on boundary conditions and structure configurations. Moreover, if we use the mode shape functions as the basis function of displacement, the mass matrix is also diagonal, which is well known as the orthogonality principle of mode shapes. Thus, it is reasonable to indicate that the proposed basis function space is an alternative method to determine the mode shape functions. The basis functions and the mode shapes are obtained in the same way: both of them are uniquely determined from boundary conditions and configurations (Jacobian matrix determinant). The basis functions and the mode shapes both have the same property: orthogonality and diagonal mass matrix. That is the reason it is called as Consistent Orthogonal Basis Functions. The basis functions are obtained consistently with mode shapes and the property of the basis functions are also consistent with mode shapes.

5

ARTICLE IN PRESS

JID: MS T. Li, R.K. Kapania

[m5GeSdc;August 18, 2017;9:38] International Journal of Mechanical Sciences 000 (2017) 1–19

Fig. 4. Orthogonal basis functions for fixed-fixed non-uniform beam.

Fig. 5. Orthogonal basis functions for fixed-free non-uniform beam.

The displacement function of a plate structure can be understood as the tensor product of the displacement functions for a 1-D beam. For a plate structure with all-clamped boundary condition, the basis functions for x, y dimension will be the tensor product of the plots shown in Fig. 2.

The Green strain tensor and its variation have the formulation as follows:

5. Nonlinear finite element implementation details

where 𝑡 𝜀𝑘𝑖𝑗 , is the Green strain tensor increment,𝑡 𝜀𝑘𝑖𝑗 =𝑡𝑡+𝑑𝑡 𝜀𝑘𝑖𝑗 −𝑡𝑡 𝜀𝑘𝑖𝑗 =𝑡𝑡+𝑑𝑡

𝑡+𝑑𝑡 𝑘 𝜀𝑖𝑗 =𝑡 𝜀𝑘𝑖𝑗 =𝑡 𝑒𝑘𝑖𝑗 +𝑡 𝜂𝑖𝑗𝑘 𝑡

𝛿𝑡𝑡+𝑑𝑡 𝜀𝑘𝑖𝑗 = 𝛿𝑡 𝜀𝑘𝑖𝑗 = 𝛿𝑡 𝑒𝑘𝑖𝑗 + 𝛿𝑡 𝜂𝑖𝑗𝑘

𝜀𝑘𝑖𝑗 ; 𝑡 𝑒𝑘𝑖𝑗 = 12 (𝑢𝑘𝑖,𝑗 + 𝑢𝑘𝑗,𝑖 ) is the linear part of 𝑡 𝜀𝑘𝑖𝑗 ; 𝑡 𝜂𝑖𝑗 = 12 𝑢𝑘𝑠,𝑖 𝑢𝑘𝑠,𝑗 is the nonlinear part of 𝑡 𝜀𝑘𝑖𝑗 . The left-low corner t in 𝑡 𝜀𝑘𝑖𝑗 ,𝑡 𝑒𝑘𝑖𝑗 and 𝑡 𝜂𝑖𝑗𝑘 means that the derivatives are calculated with respect to the time t configuration coordinate. The Kirchhoff stress tensor has the following formulation:

The displacement and the traction field for the kth layer are rewritten in matrix form as: 𝑢𝑘 = 𝑁𝑢 𝑢̃ 𝑘 𝑇 𝑘 = 𝑁𝑇 𝑇̃ 𝑘

(5 – 1a,b)

𝑡+𝑑𝑡 𝑘 𝑆𝑖𝑗 𝑡

𝑢̃ 𝑘 , 𝑇̃ 𝑘

where are vectors of DOFs (Amns , Bmns , Cmns , amn , bmn , cmn ). Nu , NT are matrices of basis functions. It is noted that for the displacement variables, the displacement increment and the displacement variation are:

=𝑡 𝑆𝑖𝑗𝑘 +𝑡𝑡 𝑆𝑖𝑗𝑘 =𝑡 𝑆𝑖𝑗𝑘 +𝑡 𝜎𝑖𝑗𝑘

(5–4)

where 𝑡 𝑆𝑖𝑗𝑘 is the increment of the Kirchhoff stress tensor based on the configuration at time t and 𝑡 𝜎𝑖𝑗𝑘 is the Cauchy stress tensor at time t. Based on the relation between Kirchhoff and the Cauchy stress tensors [1], we have: ( ) ( ) 𝑘 𝑘 𝑡 𝑘 𝑘 𝑡 𝑘 𝑡 𝑘 (5–5) 𝑡 𝑆𝑠𝑙 = 𝜎𝑠𝑙 − 𝜎𝑠𝑝 𝑡 𝑢𝑙,𝑝 − 𝜎𝑙𝑝 𝑢𝑠,𝑝

𝑢𝑘𝑖 =𝑡+𝑑𝑡 𝑢𝑘𝑖 −𝑡 𝑢𝑘𝑖 𝛿𝑢𝑘𝑖 = 𝛿 𝑡+𝑑𝑡 𝑢𝑘𝑖

(5 – 3a,b)

(5–2a,b)

Based on the Jaumann constitutive theory [1], we have (𝑡 𝑘 ) 𝑡 𝑘 (𝑡 𝑘 ) 𝑡 𝑘 (𝑡 𝑘 ) 𝑘 𝐽 𝑡 𝑆𝑖𝑗 = 𝐷𝑖𝑗𝑠𝑙 𝑒𝑠𝑙 − 𝜎𝑖𝑝 𝑒𝑗𝑝 − 𝜎𝑗𝑝 𝑒𝑖𝑝

where 𝑡+𝑑𝑡 𝑢𝑘𝑖 ,𝑡 𝑢𝑘𝑖 are the ith displacement component for the kth layer at time t + dt and time t configurations. 6

(5–6)

ARTICLE IN PRESS

JID: MS

[m5GeSdc;August 18, 2017;9:38]

T. Li, R.K. Kapania

International Journal of Mechanical Sciences 000 (2017) 1–19

(

𝐽 is the Jaumann elastic-plastic constitutive tensor and 𝑡 𝜎 𝑘 is where 𝐷𝑖𝑗𝑠𝑙 𝑖𝑗 the Cauchy stress tensor at time t. The Jaumann stress tensor is used in the constitutive model, which is necessary for large deformation analysis. The Jaumann stress tensor and its corresponding strain rate tensor are objective. It is known that the engineering strain is not applicable for large deformation analysis. A finite rigid rotation will also generate non-zero values of engineering strain, which violates the basic principle that rigid motion has no strain. An objective strain tensor means that the strain is always zero for a finite rotation motion. The Green strain tensor is objective. Taking (5–3a, b) and (5–4) into the variation of strain energy and consider (5–5, 5–6), the variation of the strain energy will be: 𝑡+𝑑𝑡

∫Ω𝑡,ℎ 𝑡

𝑆𝑖𝑗𝑘 𝛿𝑡𝑡+𝑑𝑡 𝜀𝑘𝑖𝑗 𝑑 𝑡 𝑉 =

(

∫Ω𝑡,ℎ

)( 𝑘

𝐽 𝐷𝑖𝑗𝑠𝑙 𝑡 𝑒𝑠𝑙

∫Ω𝑡,ℎ

+

∫Ω𝑡,ℎ

𝑡 𝑘 𝜎𝑖𝑗 𝛿 𝑡 𝜂𝑖𝑗𝑘

∫Ω𝑡,𝑘

−

𝑘 𝑡 𝑒𝑝𝑖

)(

𝑡+𝑑𝑡

(5–7)

A vector is defined for convenience. [ ] 𝛿⃗𝑠 = 𝛿𝑠1 , 𝛿𝑠2 , 𝛿𝑠3

(5–8) 𝐾𝑛1 =

𝐾𝑛2 =

It is noted that 𝛿⃗𝑖 𝑁𝑢 is essentially the i row of Nu . Thus, the 𝑡 𝑒𝑘𝑖𝑗 can be represented as: (( ( ) ) ( ) ) 1 1 𝑘 𝑘 𝑘 ⃗ ⃗ 𝑒 = 𝑢 + 𝑢 𝑁 + 𝑁 𝑢̃ 𝑘 = 𝛿 𝛿 𝑡 𝑖𝑗 2 𝑡 𝑖,𝑗 𝑡 𝑗,𝑖 2 𝑡 𝑖 𝑢 ,𝑗 𝑡 𝑗 𝑢 ,𝑖

(

(5–11) 𝑄=

( 𝑘 )( 𝑘 ) 𝑡 1 𝐽 𝐷𝑖𝑗𝑠𝑙 𝛿𝑡 𝑒𝑖𝑗 𝑑 𝑉 = 𝛿 𝑢̃ 𝑘𝑇 𝑡 𝑒𝑠𝑙 4 ] (( ) ( ) )𝑇 ( ( ) ( ) ) 𝐽 𝑡 ⃗ ⃗ ⃗ ⃗ 𝐷𝑖𝑗𝑠𝑙 𝑡 𝛿𝑖 𝑁𝑢 +𝑡 𝛿𝑗 𝑁𝑢 +𝑡 𝛿𝑙 𝑁𝑢 𝑑 𝑉 𝑢̃ 𝑘 𝑡 𝛿𝑘 𝑁𝑢 ,𝑖

,𝑙

,𝑠

( ) 𝛿𝑡 𝜂𝑖𝑗𝑘 𝑡 𝜎𝑖𝑗𝑘 𝑑𝑉 =

)( )]( ) [( 1 𝑡 𝑘 𝛿 𝑡 𝑢𝑘𝑠,𝑖 𝑡 𝑢𝑘𝑠,𝑗 𝜎𝑖𝑗 𝑑𝑉 ∫Ω𝑡,𝑘 2 ∫Ω𝑡,𝑘 [( )( ) ( )( )]( ) 1 𝑡 𝑘 𝛿𝑡 𝑢𝑘𝑠,𝑖 𝑡 𝑢𝑘𝑠,𝑗 + 𝑡 𝑢𝑘𝑠,𝑖 𝛿𝑡 𝑢𝑘𝑠,𝑗 = 𝜎𝑖𝑗 𝑑𝑉 ∫Ω𝑡,𝑘 2 ( [ ( ) [ ] ( ) ( )] 1 𝑡 𝑘 𝑘 ⃗ ⃗ 𝛿 𝑢̃ 𝑘 = 𝜎𝑖𝑗 ( 𝑁 ) ( 𝛿 𝑢 ̃ ) 𝑁 𝛿 𝑡 𝑠 𝑢 ,𝑖 𝑡 𝑠 𝑢 ∫Ω𝑡,𝑘 2 ,𝑗 [( ) ( )][ ( ) ( ) ]) ⃗ + 𝑡 𝛿⃗𝑠 𝑁𝑢 𝑑𝑉 𝑢̃ 𝑘 𝛿 𝑢̃ 𝑘 𝑡 𝛿𝑠 𝑁𝑢 = 𝛿 𝑢̃ 𝑘𝑇

∫Ω𝑡,𝑘

𝑡

𝛿⃗𝑠 𝑁𝑢

,𝑖

,𝑗

𝑡 𝜎𝑖𝑗

[ ) ( 𝑡

𝛿⃗𝑠 𝑁𝑢

∫Ω𝑡.𝑘

[(

𝑘 𝑡 𝑒𝑝𝑖

)(

𝑘 𝑡 𝑒𝑝𝑗

) ] ,𝑗

𝑑

) 𝑑𝑉

𝑢̃ 𝑘

)] 𝑑𝑉 =

(5–18)

,𝑖

[( 𝑡

𝛿⃗𝑖 𝑁𝑢

) ,𝑝

( ) ]𝑇 ( ) [ ( ) ( ) ] 𝑡 𝑘 +𝑡 𝛿⃗𝑝 𝑁𝑢 𝜎𝑖𝑗 𝑡 𝛿⃗𝑗 𝑁𝑢 +𝑡 𝛿⃗𝑝 𝑁𝑢 𝑑𝑡𝑉 ,𝑖

,𝑝

,𝑗

𝑡+𝑑𝑡 𝑘 𝑡 ∫𝑡+𝑑𝑡 Γ± 𝑡+𝑑𝑡 𝑠

(

𝛿⃗𝑠 𝑁𝑢

)𝑇

𝑑 𝑡+𝑑𝑡 𝑆

(

) +

∫Ω𝑡+𝑑𝑡,𝑘

𝑡+𝑑𝑡 𝑘 𝑓 𝑡+𝑑𝑡 𝑠

(5–20) (

𝛿⃗𝑠 𝑁𝑢

)𝑇

) 𝑑 𝑡+𝑑𝑡 𝑉

𝑡+𝑑𝑡

⎡ ⃗𝑖 ⎢ 𝜕 𝑡+𝑑𝑡 𝑥 (𝑡+𝑑𝑡 ) 𝑡+𝑑𝑡 1 𝑆⃗ = 𝑡+𝑑𝑡 𝑛⃗ 𝑑 𝑆 = det ⎢ 𝜕 𝜉1 ⎢ 𝜕 𝑡+𝑑𝑡 𝑥 1 ⎢ ⎣ 𝜕 𝜉2

𝑗⃗

𝜕 𝑡+𝑑𝑡 𝑥2 𝜕 𝜉1 𝜕 𝑡+𝑑𝑡 𝑥2 𝜕 𝜉2

+ 𝑢𝑘𝑖

⃗ 𝑘

⎤ ⎥

𝜕 𝑡+𝑑𝑡 𝑥3 ⎥𝑑 𝜉 𝑑 𝜉 𝜕 𝜉1 ⎥ 1 2 𝜕 𝑡+𝑑𝑡 𝑥3 ⎥ 𝜕 𝜉2 ⎦

(5–23)

(5–24)

⃗ are the basis vector of space rectangular coordinate system; where ⃗𝑖, 𝑗⃗, 𝑘 𝜉 i ,i = 1, 2, 3 are the natural coordinates; t + dt xi is the configuration i coordinate at time t + dt; t xi is the configuration i coordinate at time t. Thus, the area nominal vector is:

( )( ) 2𝑡 𝜎𝑖𝑗𝑘 𝛿𝑡 𝑒𝑘𝑝𝑖 𝑡 𝑒𝑘𝑝𝑗 𝑑𝑉

∫Ω𝑡,𝑘 [( ) ( ) ]𝑇 ( ) 1 𝑡 𝑘 = 𝛿 𝑢̃ 𝑘𝑇 𝜎𝑖𝑗 𝛿⃗𝑖 𝑁𝑢 +𝑡 𝛿⃗𝑝 𝑁𝑢 𝑡 ∫Ω𝑡,𝑘 2 ,𝑝 ,𝑖 ] ) [( ) ⃗ +𝑡 (𝛿⃗𝑝 𝑁𝑢 ),𝑗 𝑑𝑉 𝑢̃ 𝑘 𝑡 𝛿𝑗 𝑁𝑢

1 2 ∫Ω𝑡,𝑘

𝑡+𝑑𝑡 𝑘 𝑡 𝑘 𝑥𝑖 = 𝑥𝑖

(5–13)

⎡ 𝜕 𝑡+𝑑𝑡 𝑥𝑘2 ⎢ 𝜕 𝜉1 ⎢ 𝑡+𝑑𝑡 𝑘 𝜕 𝑥3 𝑡+𝑑𝑡 ⃗ 𝑆 =⎢ 𝑑 ⎢ 𝜕 𝜉1 ⎢ 𝜕 𝑡+𝑑𝑡 𝑥𝑘 1 ⎢ ⎣ 𝜕 𝜉1

(

,𝑝

∫Ω𝑡,𝑘

[( )𝑇 ]( )[ ] 𝑡 ⃗ 𝜎𝑖𝑗 𝑡 (𝛿⃗𝑠 𝑁𝑢 ),𝑗 𝑑𝑉 𝑡 𝛿𝑠 𝑁 𝑢

)

where p is the fluid pressure; ns is the sth component of the unit normal vector pointing outside 𝑛⃗, 𝑛⃗ = [𝑛1 , 𝑛2 , 𝑛3 ]𝑇 ; t + dt refers to the current time. Then, assuming the fluid pressure is applied on the surface 𝜉 3 = 1

The third term in the strain energy variation is: 𝑡 𝑘 𝜎𝑖𝑗 𝛿

(5–16)

The load tangential stiffness matrix is also considered, for example, for the case of fluid pressure that is always normal to the deformed structure’s outer surface: )𝑇 (𝑡+𝑑𝑡 )( ( )𝑇 (𝑡+𝑑𝑡 ) 𝑡+𝑑𝑡 𝑁𝑢 𝑄= 𝑛 𝛿⃗𝑠 𝑁𝑢 𝑝𝑑 𝑡+𝑑𝑡 𝑆 = 𝑛 𝑝𝑑 𝑆 𝑡+𝑑𝑡 ∫𝑡+𝑑𝑡 Γ𝜎 𝑡+𝑑𝑡 𝑠 ∫𝑡+𝑑𝑡 Γ𝜎 (5–22)

The second term in the strain energy variation is:

)𝑇 ](

,𝑖

(5–21)

(5–12)

,𝑖

,𝑗

) ( ) )]𝑇 +𝑡 𝛿⃗𝑗 𝑁𝑢 𝑑𝑉

And the internal/external force vector are: ( ) [( )( ( ) ( ) )]𝑇 1 𝑡 𝑘 ⃗ ⃗ 𝐹 = 𝜎𝑖𝑗 𝑡 𝛿𝑖 𝑁𝑢 +𝑡 𝛿𝑗 𝑁𝑢 𝑑𝑉 ∫Ω𝑡,𝑘 2 ,𝑗 ,𝑖

(5–10)

Thus, the first term in the variation of the strain energy is:

[(

)

(5–19)

Similarly, the variation of the linear strain tensor is: (( ( ) ) ( ) ) 1 1 ⃗ ⃗ 𝛿𝑡 𝑒𝑘𝑖𝑗 = 𝛿 𝑡 𝑢𝑘𝑖,𝑗 +𝑡 𝑢𝑘𝑗,𝑖 = 𝑁 + 𝑁 𝛿 𝑢̃ 𝑘 𝛿 𝛿 2 2 𝑡 𝑖 𝑢 ,𝑗 𝑡 𝑗 𝑢 ,𝑖

(

𝑡

𝛿⃗𝑖 𝑁𝑢

𝑡+𝑑𝑡 𝑘 𝑘 𝑡+𝑑𝑡 𝑡+𝑑𝑡 𝑘 𝑘 𝑡+𝑑𝑡 𝑡 𝛿𝑢 𝑑 𝑆+ 𝑓 𝛿𝑢 𝑑 𝑉 ∫𝑡+𝑑𝑡Γ+ 𝑡+𝑑𝑡 𝑠 𝑠 ∫Ω𝑡+𝑑𝑡,𝑘 𝑡+𝑑𝑡 𝑠 𝑠 ( ) ( )𝑇 𝑡+𝑑𝑡 𝑘 ⃗ 𝑡+𝑑𝑡 = 𝛿 𝑢̃ 𝑘𝑇 𝑡 𝑁 𝑑 𝑆 𝛿 𝑠 𝑢 𝑠 ∫𝑡+𝑑𝑡 Γ+ 𝑡+𝑑𝑡 ( ) ( )𝑇 𝑡+𝑑𝑡 𝑘 ⃗ 𝑡+𝑑𝑡 +𝛿 𝑢̃ 𝑘𝑇 𝑓 𝑁 𝑑 𝑉 𝛿 ∫Ω𝑡+𝑑𝑡,𝑘 𝑡+𝑑𝑡 𝑠 𝑠 𝑢

(5–9)

,𝑗

)( (

(5–17)

𝑢𝑘𝑠 = 𝛿⃗𝑠 𝑁𝑢 𝑢̃ 𝑘 , (𝑠 = 1, 2, 3)

∫𝑡 𝑉

𝑡 𝑘 𝜎𝑖𝑗

𝑊𝑘 =

(

Thus, we have:

∫𝑡 𝑉 [

[(

where 𝑓𝑠𝑘 is the body force and 𝑡𝑘𝑠 is the external traction force for the kth layer. Then the tangential stiffness matrices are obtained: [ ] (( ) ( ) )𝑇 ( ( ) ( ) ) 1 𝐽 𝑡 ⃗ ⃗ ⃗ ⃗ 𝐾𝑙 = 𝐷𝑖𝑗𝑠𝑙 𝑁 + 𝑁 𝑁 + 𝑁 𝑑 𝑉 𝛿 𝛿 𝛿 𝛿 𝑡 𝑖 𝑢 𝑡 𝑗 𝑢 𝑡 𝑘 𝑢 𝑡 𝑙 𝑢 4 ∫𝑡 𝑉 ,𝑗 ,𝑖 ,𝑙 ,𝑠

)] 𝑘 𝑑𝑡𝑉 𝑡 𝑒𝑝𝑗

𝑡 𝑘 𝜎𝑖𝑗 𝛿𝑡 𝑒𝑘𝑖𝑗 𝑑 𝑡 𝑉

∫Ω𝑡,𝑘

1 2

For the body force term and external distributed force, we have:

𝛿𝑡 𝑒𝑘𝑖𝑗 𝑑 𝑡 𝑉 (

= 𝛿 𝑢̃ 𝑘𝑇

(5–15)

)

[

+

𝑡 𝑘 𝜎𝑖𝑗 𝛿𝑡 𝑒𝑘𝑖𝑗 𝑑𝑉

(5–14)

7

𝜕 𝑡+𝑑𝑡 𝑥𝑘3 𝜕 𝜉2

𝜕 𝑡+𝑑𝑡 𝑥𝑘1 𝜕 𝜉2 𝜕 𝑡+𝑑𝑡 𝑥𝑘2 𝜕 𝜉2

− − −

𝜕 𝑡+𝑑𝑡 𝑥𝑘2 𝜕 𝑡+𝑑𝑡 𝑥𝑘3 ⎤ 𝜕 𝜉2 𝜕 𝜉1 ⎥ ⎥ 𝜕 𝑡+𝑑𝑡 𝑥𝑘1 𝜕 𝑡+𝑑𝑡 𝑥𝑘3 ⎥ 𝜕 𝜉1 𝜕 𝑡+𝑑𝑡 𝑥𝑘1 𝜕 𝜉2

𝜕 𝜉2 𝜕 𝑡+𝑑𝑡 𝑥𝑘2 𝜕 𝜉1

⎥ ⎥ ⎥ ⎦

𝑑 𝜉1 𝑑 𝜉2

ARTICLE IN PRESS

JID: MS T. Li, R.K. Kapania (

International Journal of Mechanical Sciences 000 (2017) 1–19 ) (

)

⎡ 𝜕 𝑡 𝑥𝑘2 +𝑢𝑘2 𝜕 𝑡 𝑥𝑘3 +𝑢𝑘3 − ⎢ 𝜕𝜉 ⎢ ( 1 ) ( 𝜕 𝜉2 ) ⎢ 𝜕 𝑡 𝑥𝑘3 +𝑢𝑘3 𝜕 𝑡 𝑥𝑘1 +𝑢𝑘1 =⎢ − ⎢ ( 𝜕 𝜉1 ) ( 𝜕 𝜉2 ) ⎢ 𝜕 𝑡 𝑥𝑘 +𝑢𝑘 𝜕 𝑡 𝑥𝑘 +𝑢𝑘 1 1 2 2 ⎢ − 𝜕 𝜉2 ⎣ 𝜕 𝜉1

( ) ( ) 𝜕 𝑡 𝑥𝑘2 +𝑢𝑘2 𝜕 𝑡 𝑥𝑘3 +𝑢𝑘3 ⎤

⎥ ( ) ( )⎥ 𝜕 𝑡 𝑥𝑘1 +𝑢𝑘1 𝜕 𝑡 𝑥𝑘3 +𝑢𝑘3 ⎥ ⎥𝑑 𝜉1 𝑑 𝜉2 𝜕 𝜉1 𝜕 𝜉2 ( ) ( )⎥ 𝜕 𝑡 𝑥𝑘1 +𝑢𝑘1 𝜕 𝑡 𝑥𝑘2 +𝑢𝑘2 ⎥ ⎥ 𝜕 𝜉2 𝜕 𝜉1 ⎦ 𝜕 𝜉2

Table 1a The comparison of the axial displacement and the shear stress.

𝜕 𝜉1

(5–25)

𝑢̄ = 𝐿𝑡 3 𝑢( 𝐿4 , 𝑡) 𝜏̄ = 𝜏𝑥𝑧 ( 𝐿4 , 𝑧1 )∕𝑝0 2

1

1

∫−1 ∫−1

(5–26a)

𝑢̄ = 𝜏̄ =

1

1

𝑝𝑁𝑢𝑇 𝑁𝑛 𝑑 𝜉2 𝑑 𝜉3

(5-28a)

where 𝛿⃗𝑖 is defined in (5–8, 5–9) and 𝑘

⎡𝑇1 ⎤ 𝑇 = ⎢𝑇2 ⎥ ⎢ ⎥ ⎣𝑇3 ⎦

(5-28b)

Taking the interpolation matrix of traction field defined in (5–1b) into (5–28a): 𝑇 𝑘 = 𝑁𝑇 𝑇̃ 𝑘

(5–29)

where 𝑇̃ is the degree-of-freedom vector for traction field. Thus, consider (5–28), we have: ( ) 𝑘 𝑇𝑖 = 𝛿⃗𝑖 𝑁𝑇 𝑇̃ 𝑘

(5–30)

2

𝑡,𝑘

𝑛𝑙

𝑘=1 𝑛𝑙

𝑘=1

5.623 × 10 −2.025

−7

Proposed method 5.329 × 10 − 7 −2.077

transverse shear stress are defined as: 𝑢̄ = 𝐿𝑡 3 𝑢(𝑥, 𝑧) and 𝜏̄ = 𝜏𝑥𝑧 (𝑥, 𝑧)∕𝑝0 , where p0 = 1 × 106 Pa. The results are compared in the Table 1a. In Table 1a, the nondimensional axial displacement and the shear stress are presented. The proposed method has a matching results with 3-D finite element solution and the RTZ solution presented in [38]. The plot of the non-dimensional shear stress along the thickness direction is given in Fig. 6. In order to have the distribution of the transverse shear stress along the thickness direction, we use multiple finite element layers to represent one material layer. As Fig. 6 shows, each point represents an interface. Problem 2: The problem 2 is also taken from the section 5.2 of [38]. For this case, the loads are same as those of Problem 1. The thicknesses for the 5 layers are [0.025/0.025/0.001 /0.025/0.025]. The Young’s modulus are [73/21.9/0.073/73/21.9]GPa and the shear modulus are [29.2/8.76/0.029/29.2/8.76]. Other parameters and boundary conditions are same as those in Problem 1. The results are compared in Table 1b. It is noted that the results also match well.

Thus, for the last two terms of the weak form (3–1), the discrete formulations will be: 𝑛𝑙 ⟨ ⟩ ∑ 𝛿𝑢𝑘𝑖 , 𝑇𝑖𝑘 Γ± 𝑡,𝑘 𝑘=1 ( ( )𝑇 𝑘 ( ) ) 𝑛𝑙 ∑ 𝛿⃗𝑖 𝑁𝑇 𝑑𝑆 𝑇̃ 𝑘 = 𝛿 𝑢̃ 𝑘𝑖 ∫Γ± 𝛿⃗𝑖 𝑁𝑢 (5–31) 𝑡,𝑘 𝑘=1 ( ) 𝑛𝑙 ( )𝑇 ( ) ∑ = 𝛿 𝑢̃ 𝑘𝑖 ∫Γ± 𝑁𝑢 𝑁𝑇 𝑑𝑆 𝑇̃ 𝑘 ∑ ⟨ 𝑘 𝑘⟩ 𝛿𝑇𝑖 , 𝑢𝑖 Γ± ( 𝑡,𝑘( )𝑇 ( ) ) ∑ ̃𝑘 = 𝛿 𝑇 ∫Γ± 𝛿⃗𝑖 𝑁𝑇 𝛿⃗𝑖 𝑁𝑢 𝑑𝑆 𝑢̃ 𝑘𝑖 𝑘=1 ( 𝑡,𝑘 ) 𝑛𝑙 ( )𝑇 ( ) ∑ 𝑘 ̃ = 𝛿 𝑇 ∫Γ± 𝑁𝑇 𝑁𝑢 𝑑𝑆 𝑢̃ 𝑘𝑖

5.867 × 10 −2.023

3-D finite element −7

In this section, the proposed method is used to solve several problems that have been analyzed by other well-established methods in the previously publicized papers. We will compare the results to verify the accuracy of this paper’s method, especially in predicting the interface shear stress field. Problem 1: The problem 1 is taken from the Section 5.2 of [38], in which the clamped-clamped laminated beam under sinusoidal pressure and shear traction is analyzed. The laminated beam is also a branch of laminated plate, in which the y dimension is fixed. The length of the beam is L = 1m. The plane stress assumption is applied. The thicknesses for the 3 layers are [0.01/0.08/0.01] from the bottom layer to the top layer. The top layer is applied by a normal pressure of 𝑃̂𝑡 = −0.5 sin(𝜋𝑥) MPa and a shear traction of 𝑇̂𝑡 = −4 cos(𝜋𝑥) MPa; the bottom layer is applied by a normal pressure of 𝑃̂𝑡 = 0.5 sin(𝜋𝑥) MPa and a shear traction of 𝑇̂𝑡 = 2 cos(𝜋𝑥) MPa. The Young’s modulus for the 3 layers are [73/0.073/21.9]GPa and the shear modulus for the 3 layers are [29.2/0.029/8.76]GPa. The isotropic elastic constitutive theory is applied. The non-dimensional axial displacement and the non-dimensional

(5–27b)

𝑘

𝑢( 𝐿4 , 𝑡) 𝜏𝑥𝑧 ( 𝐿4 , 𝑧1 )∕𝑝0

6.1. Verification

(5–27a)

where p is the distributed pressure function. For the ith component of the traction field, we have: 𝑇𝑖𝑘 = 𝛿⃗𝑖 𝑇 𝑘

𝑡2 𝐿3

The numerical tests results are presented in this section. In Section 6.1, a verification study is presented. In Section 6.2, the sandwich plate under distributed load is studied. In Section 6.3, the sandwich plate under concentrated force is analyzed. The interface stress and the displacement are both solved and compared with ANSYS.

[ 𝑡 𝑘( ) ( ) ] [ 𝜕 𝑡 𝑥𝑘 ( ) ) ] 𝑡 𝑘( 𝜕𝑡 𝑥𝑘3 ⎡ 𝜕 𝑥2 𝛿⃗ 𝑁 ⎤ 2 ⃗2 𝑁𝑢 ⃗3 𝑁𝑢 + 𝜕 𝑥3 𝛿⃗2 𝑁𝑢 + − 𝛿 𝛿 3 𝑢 𝜕 𝜉3 𝜕 𝜉3 𝜕 𝜉2 ⎢ 𝜕 𝜉2 ,3 ,2 ,2 ,3 ⎥ [ ] [ ] ⎢ 𝜕 𝑡 𝑥𝑘 ( ) ( ) ( ) ⎥ 𝜕 𝑡 𝑥𝑘 𝜕 𝑡 𝑥𝑘 𝜕 𝑡 𝑥𝑘 ⎥ 𝑁𝑛 = ⎢ 𝜕𝜉 3 𝛿⃗1 𝑁𝑢 + 𝜕 𝜉 1 𝛿⃗3 𝑁𝑢 − 𝜕 𝜉 1 (𝛿⃗3 𝑁𝑢 ),3 + 𝜕 𝜉 3 𝛿⃗1 𝑁𝑢 ⎢ 2 3 2 3 ,3 ,2 ,2 ⎥ ] [ 𝑡 𝑘( ]⎥ ) ) ( ) ⎢ [ 𝜕 𝑡 𝑥𝑘 ( 𝑘 𝑘 𝑡 𝑡 𝜕𝑥 𝜕𝑥 𝜕𝑥 1 ⃗ 𝑢 ⎢ ⎥ + 2 (𝛿⃗ 𝑁 ) − 𝜕 𝜉 2 𝛿⃗1 𝑁𝑢 + 𝜕 𝜉 1 𝛿𝑁 𝛿⃗ 𝑁 ⎣ 𝜕 𝜉2 2 𝑢 , 3 𝜕 𝜉3 1 𝑢 , 2 2 3 ,3 ,2 ⎦

𝑘=1

Proposed method 2.787 × 10 − 6 −2.372

6. Numerical tests

If the water pressure is applied on the surface 𝜉 2 − 𝜉 3 , the load tangential stiffness matrix will be: ∫−1 ∫−1

3-D finite element 2.756 × 10 − 6 −2.316

RZT

𝑝𝑁𝑢𝑇 𝑁𝑛 𝑑 𝜉1 𝑑 𝜉2

[ ( ) ( ) ] [ 𝜕 𝑡 𝑥𝑘 ( ) ( ) ] 𝜕 𝑡 𝑥𝑘3 𝜕 𝑡 𝑥𝑘3 ⎡ 𝜕 𝑡 𝑥𝑘2 ⃗ ⎤ 2 ⃗ ⃗ ⃗ ⎢ 𝜕𝜉1 𝛿3 𝑁𝑢 ,2 + 𝜕 𝜉2 𝛿2 𝑁𝑢 ,1 − 𝜕 𝜉2 𝛿3 𝑁𝑢 ,1 + 𝜕 𝜉1 𝛿2 𝑁𝑢 ,2 ⎥ ⎢[ ⎥ ) ( ) ] [ 𝜕 𝑡 𝑥𝑘 ( ) ( ) ]⎥ 𝜕 𝑡 𝑥𝑘 𝜕 𝑡 𝑥𝑘 ⎢ 𝜕 𝑡 𝑥𝑘 ( 𝑁𝑛 = ⎢ 𝜕 𝜉 3 𝛿⃗1 𝑁𝑢 + 𝜕 𝜉 1 𝛿⃗3 𝑁𝑢 − 𝜕 𝜉 1 𝛿⃗3 𝑁𝑢 + 𝜕 𝜉 3 𝛿⃗1 𝑁𝑢 ⎥ 2 1 2 ,2 ,1 ,2 ,1 ⎢[ 1 ] [ 𝑡 𝑘( ]⎥⎥ ) ( ) ) ( ) 𝑘 𝑘 𝑡 𝑡 ⎢ 𝜕 𝑡 𝑥𝑘 ( 𝜕𝑥 𝜕𝑥 𝜕𝑥 ⎢ 𝜕 𝜉 1 𝛿⃗2 𝑁𝑢 + 𝜕 𝜉 2 𝛿⃗1 𝑁𝑢 ⎥ − 𝜕 𝜉 2 𝛿⃗1 𝑁𝑢 + 𝜕 𝜉 1 𝛿⃗2 𝑁𝑢 1 2 1 2 ⎣ ,2 ,1 ,2 ,1 ⎦ (5–26b)

𝐾𝑄 =

RZT 2.763 × 10 − 6 −2.323

Table 1b The comparison of the axial displacement and the shear stress.

The linear part from (5–22) can be picked up as the load tangential stiffness matrix, denoted as KQ , thus 𝐾𝑄 =

[m5GeSdc;August 18, 2017;9:38]

(5–32)

𝑡,𝑘

8

ARTICLE IN PRESS

JID: MS T. Li, R.K. Kapania

[m5GeSdc;August 18, 2017;9:38] International Journal of Mechanical Sciences 000 (2017) 1–19

Table 1c The comparison of the max deflection and the interface shear stresses. a/t

Method

𝑤̄

𝜏̄𝑦𝑧

𝜏̄𝑥𝑧

4

Present Shi [39] ABAQUS Pagano [40,41] Di Sciuva [42] Lee [43] Mantari [44]

2.0057 1.9734 2.0061 2.0042 1.8910 1.9538 1.9434

0.2247 0.2297 0.2226 0.2172 0.2542 0.2336 0.2010

0.2651 0.2409 0.2560 0.2560 0.2620 0.2390 0.2450

10

Present Shi [39] ABAQUS Pagano [40,41] Di Sciuva [42] Lee [43] Mantari [44]

0.7537 0.7549 0.7522 0.7528 0.7233 0.7529 0.7342

0.1272 0.1239 0.1377 0.1228 0.1230 0.1243 0.1150

0.3698 0.3566 0.3581 0.3570 0.3643 0.3561 0.3140

20

Present Shi [39] ABAQUS Pagano [40,41] Di Sciuva [42] Lee [43] Mantari [44]

0.5172 0.5176 0.5172 0.5166 0.5076 0.5165 0.5113

0.0974 0.0940 0.0966 0.0938 0.0933 0.0940 0.0900

0.3981 0.3845 0.3848 0.3850 0.3869 0.3845 0.3310

100

Present Shi [39] ABAQUS Pagano [40,41] Di Sciuva [42] Lee [43] Mantari [44]

0.4356 0.4356 0.4362 0.4351 0.4343 0.4347 0.4353

0.0864 0.0829 0.0813 0.0828 0.0828 0.0828 0.0810

0.4093 0.3946 0.3935 0.3950 0.3948 0.3947 0.3370

Fig. 6. The distribution of the shear stress 𝜏̄𝑥𝑧 along the thickness direction.

6.2. Distributed load A sandwich plate structure with 3 different boundary conditions is analyzed in this Section, as Fig. 8 shows. Throughout this paper, we assume that the elastic material is used. The Young’s Modulus is defined as Es = 206GPa for the top and bottom layer and Ec = 2.06GPa for the core layer. The Possion ratio is given as 𝜈 = 0.3. The thicknesses of the top and the bottom layers are given as ts = 0.01m and the core layer is defined as tc = 0.1m. The sandwich plate’s length and width are given as a = b = 1m. The water pressure is applied on the top surface of the sandwich plate. The direction of the water pressure is a function of the deformed sandwich plate, thus causing a nonlinear problem due to loads. The geometric nonlinearity is also considered because of large deformation analysis. All the numerical results are verified by using a traditional finite element software: ANSYS. The second order Lagrangian element Solid 186 element is used. The element initial shape is a cube. The meshing density will be discussed in Section 6.2.1–6.2.3, where the convergence study for meshing density is presented. On the interfaces between the core layer and top/bottom layers, a conforming mesh is used such that the displacement is coupled. They share the same nodes on the interface. The displacement of the core layer and the adjacent top/bottom layers is assumed to be the same. The delamination will be considered in future work. The displacement and the interface shear stress T1 ,T2 are analyzed by using the proposed method. Results are compared with those obtained from ANSYS.

Fig. 7. The distribution of the shear stress along the thickness direction.

The non-dimensional shear stress is plotted against the thickness coordinate in Fig. 7. Problem 3: The problem 3 is taken from [39], in which the laminated plate with cross-ply [0/90/0] core layer is analyzed. The 3 layers have the same ∑ thickness tk (k = 1, 2, 3) = 1 and the total thickness is 𝑡 = 3𝑖=1 𝑡𝑖 = 3. The length and width of the plate are defined as a = b = 12, 30, 60, 300. The simply supported boundary conditions are used for all the four edges. The sinusoidal pressure q = q0 sin (𝜋x/a)sin (𝜋y/b) is applied on the top surface. The orthotropic elastic constitutive model is applied here: 𝐸1 = 174.6𝐺𝑃 𝑎; 𝐸2 = 𝐸3 = 7𝐺𝑃 𝑎 𝐺12 = 𝐺13 = 3.5𝐺𝑃 𝑎; 𝐺23 = 1.4𝐺𝑃 𝑎 𝜈12 = 𝜈23 = 𝜈13 = 0.25 The

non-dimensional

𝐸 𝑡3 𝑤(𝑎∕2, 𝑏∕2, 0) 𝑞 2𝑎4 0

parameters

are

defined

as

𝑤̄ =

× 100 and

𝜏̄𝑦𝑧 = 𝜏𝑦𝑧 (𝑎∕2, 0, 0) 𝑞 𝑡 𝑎 , 𝜏̄𝑥𝑧 = 𝜏𝑥𝑧 (0, 𝑏∕2, 0) 𝑞 𝑡 𝑎 . The results of the maxi0

6.2.1. Clamped sandwich plate Since all the boundaries are clamped, we could define the displacement components:

0

mum deflection, interface shear stresses are compared in the Table 1c. In Table 1c, the ratio of length-to-thickness is defined as 4, 10, 20 and 100. Through the comparison, it is found that the proposed method shows a matching result with those previous papers [39–44] for both thin plates and thick plates.

𝜉

𝜉

𝜉

𝑢𝑘1 = 𝑅0 𝑓𝑢11 𝑓𝑢12 𝑓𝑢13 9

𝑚0 𝑛0 𝑠0 ∑ ∑∑ 𝑚=1 𝑛=1 𝑠=1

( ) 𝑛 ( ) 𝑠 ( 𝑘) 𝐴𝑚𝑛𝑠 𝜙𝑚 𝑢 𝜉1 𝜓𝑢 𝜉2 𝜑𝑢1 𝜉3 1

1

ARTICLE IN PRESS

JID: MS T. Li, R.K. Kapania

[m5GeSdc;August 18, 2017;9:38] International Journal of Mechanical Sciences 000 (2017) 1–19

Fig. 8. The boundary conditions for the laminated plates.

Fig. 9. The Consistent Orthogonal Basis Functions for 𝜉 1 .

𝜉

𝜉

𝜉

𝑢𝑘2 = 𝑅0 𝑓𝑢21 𝑓𝑢22 𝑓𝑢23 𝜉

𝜉

𝜉

𝑢𝑘3 = 𝑅0 𝑓𝑢31 𝑓𝑢32 𝑓𝑢33

𝑚0 𝑛0 𝑠0 ∑ ∑∑ 𝑚=1 𝑛=1 𝑠=1 𝑚0 𝑛0 𝑠0

∑ ∑∑

𝑚=1 𝑛=1 𝑠=1

Fig. 10. The Consistent Orthogonal Basis Functions for 𝜉 3 .

( ) 𝑛 ( ) 𝑠 ( 𝑘) 𝐵𝑚𝑛𝑠 𝜙𝑚 𝑢 𝜉1 𝜓𝑢 𝜉2 𝜑𝑢 𝜉3 2

2

2

( ) 𝑛 ( ) 𝑠 ( 𝑘) 𝐶𝑚𝑛𝑠 𝜙𝑚 𝑢 𝜉1 𝜓𝑢 𝜉2 𝜑𝑢 𝜉3 3

3

3

(6–2–1 a,b,c)

where ( ) 𝜉 𝜉 𝜉 𝑓𝑢11 = 𝑓𝑢21 = 𝑓𝑢31 = 𝜉1 1 − 𝜉1 ( ) 𝜉 𝜉 𝜉 𝑓𝑢12 = 𝑓𝑢22 = 𝑓𝑢32 = 𝜉2 1 − 𝜉2 𝜉𝑘

𝜉𝑘

𝜉𝑘

𝑓𝑢13 = 𝑓𝑢23 = 𝑓𝑢33 = 1

(6–2–2 a,b,c)

Thus, based on Section 4, the orthogonal basis function space for the sandwich plate structure will be: 𝜉 𝜑𝑢1𝑖 (𝑖 = 1, 2, 3) is the orthogonal polynomial series with weight function 2 𝜉1 (1 − 𝜉1 )2 ,

Fig. 11. The Consistent Orthogonal Basis Functions for 𝜉 1 and 𝜉 2 .

𝜉

𝜓𝑢𝑖2 (𝑖 = 1, 2, 3) is the orthogonal polynomial series with weight function 𝜉22 (1 − 𝜉2 )2 , 𝜉 𝜙𝑢1𝑖 (𝑖

Thus, the orthogonal basis function space for the sandwich plate structure will be: 𝜉 𝜑𝑢1𝑖 (𝑖 = 1, 2, 3) is the orthogonal polynomial series with weight function (𝜉 1 (1 − 𝜉 1 ))4 , 𝜉 𝜓𝑢𝑖2 (𝑖 = 1, 2, 3) is the orthogonal polynomial series with weight function (𝜉 2 (1 − 𝜉 2 ))4 , 𝜉 𝜙𝑢1𝑖 (𝑖 = 1, 2, 3) is the orthogonal polynomial series with weight function 1. The Consistent Orthogonal Basis Functions for 𝜉 1 and 𝜉 2 are plotted in Fig. 11. The difference between Figs. 9 and 11 are the assumption of shear deformation. If the plate is shear deformable, even if we apply a clamp boundary condition, the rotation at the boundaries is still non-zero. If

= 1, 2, 3) is the orthogonal polynomial series with weight function

1.

In this case, the Consistent Orthogonal Basis Functions for 𝜉 1 and 𝜉 2 are plotted in Fig. 9. For 𝜉 3 , the Consistent Orthogonal Basis Functions are plotted in Fig. 10. Another choice is also presented. The boundary conditions are defined as: ( )2 𝜉 𝜉 𝜉 𝑓𝑢11 = 𝑓𝑢31 = 𝑓𝑢31 = 𝜉12 1 − 𝜉1 ( )2 𝜉 𝜉 𝜉 𝑓𝑢11 = 𝑓𝑢21 = 𝑓𝑢31 = 𝜉12 1 − 𝜉1 𝜉

𝜉

𝜉

𝑓𝑢13 = 𝑓𝑢23 = 𝑓𝑢33 = 1

(6–2–3 a,b,c) 10

ARTICLE IN PRESS

JID: MS T. Li, R.K. Kapania

[m5GeSdc;August 18, 2017;9:38] International Journal of Mechanical Sciences 000 (2017) 1–19

Table 2d Convergence study on max deflection of mesh density for ANSYS Solid 186. DOF unit:mm

1476 21.357

2688 21.854

3429 22.758

4260 22.203

5181 22.762

Fig. 12. The deflection curves under several water pressure values. Table 2a Convergence study on max deflection of (m0 ,n0 ) for (k0 = 4). PO unit:mm

(5,5) 20.963

(6,6) 21.047

(7,7) 23.137

(8,8) 23.219

(9,9) 22.794

(10,10) 22.791

(11,11) 22.988 Fig. 13. The z-displacement (deflection) contour for the sandwich plate from ANSYS.

∗

PO is polynomial order for the displacement basis functions (a,b) refers to the displacement order in the plate x-y direction ∗ the number of DOF is calculated as ndof = 3abnl, nl = 3 ∗

Table 3a Convergence study of in-plane displacement order for shear stress (DFEM).

Table 2b Convergence study on max deflection of mesh density for ANSYS Solid 186. DOF unit:mm

9243 22.859

10,752 22.593

12,375 22.865

14,112 22.667

15,963 22.872

17,928 22.719

PO unit:MPa

45,708 22.846

(5,5) 20.980

(6,6) 21.052

(7,7) 23.138

(8,8) 23.182

(9,9) 22.772

(6,6) 65.822

(7,7) 82.202

(8,8) 71.542

(9,9) 71.174

(10,10) 70.617

(11,11) 70.062

Table 3b Convergence study of mesh density for shear stress (ANSYS Solid 186).

Table 2c Convergence study on max deflection of (m0 , n0 ) for (𝑘0 = 3). PO unit:mm

(5,5) 86.460

(10,10) 22.773

DOF unit:MPa

(11,11) 22.997

we assume the plate is shear undeformable, there will be no rotation at the boundaries. In this case, the basis functions shown in Fig. 11 will be used. The deformation history of the sandwich plate under several pressure levels is plotted in Fig. 12. The convergence study is presented in Table 2a for the displacement order of x, y dimensions. The maximum deflection is calculated under 30 MPa water pressure. In Table 2a, we assume that the thickness displacement order k0 is 4. In Table 2b, the results from ANSYS using the Solid 186 element are given. It is noted that the mesh division in the thickness direction is one for the top/bottom layer and two for the core layer. From the results in Table 2a, it is observed that, when the in-plane displacement order is greater than 7, the change of the result is less than 1% (22.988 × (1 ± 1%)). Thus, the minimum number of DOF for the proposed method will be 3 × 7 × 7 × 4 × 3 = 1764. For the traditional finite element method, from Table 2b, we can observe that the change in the result is less than 1% (22.846 × (1 ± 1%)) when the number of DOF becomes greater than 10,752. Next, we prescribe that the thickness displacement order k0 to be 3 for the proposed method and the results are given in Table 2c. To obtain the results using ANSYS, the mesh division in the thickness direction is one for the top/bottom layer and one for the core layer. The ANSYS results are given in Table 2d. From Table 2d, it is observed that the convergence condition is still not satisfied even when the number of DOF is 4260. For the proposed method, see the results in Table 2c, it is observed that the convergence

21,600 78.600

26,268 78.200

45,360 76.700

68,352 74.300

90,720 72.300

123,984 72.800

246,480 72.500

solution is obtained when the in-plane displacement order is greater than 7. The number of DOF is 3 × 7 × 7 × 3 × 3 = 1323. It is observed that, from Table 2d, the number of DOF is much lower than that in Table 2b. In Table 2b, we define the mesh size for all layers as tc /2. In Table 2d, we define the mesh size for all layers as tc . Thus, the mesh density is much larger in Table 2b. The motivation is that we prefer the element division, of which the element is close to a cube, instead of a thin cuboid. The displacement contour obtained from ANSYS and the proposed method are shown in Figs. 13 and 14, respectively. The in-plane shear stress between adjacent layers is analyzed. The convergence study for the shear stress is presented here. By increasing the number of DOF, the maximum shear stress T1 is calculated in Table 3a–b. The traction T1 is plotted against the X coordinate, for different inplane displacement order, in Fig 15. From the results, it is noted that the in-plane displacement order should be greater than 10 for an accurate solution, especially for the region near the boundaries where the shear stress solution appears to have a sharp gradient.Fig. 15(a) and (b) The shear traction T1 ’s contour is presented in Fig. 16(a) and (b). By using ANSYS, the shear traction contour is given in Fig. 16(c) and (d)). The shear traction solution is compared with that from ANSYS using a very fine mesh in Fig. 17. The normal traction is compared. It is noted that the normal traction should be 0 at the boundaries because the boundaries are all clamped. Thus, the boundary condition functions for traction field are defined 11

JID: MS

ARTICLE IN PRESS

T. Li, R.K. Kapania

[m5GeSdc;August 18, 2017;9:38] International Journal of Mechanical Sciences 000 (2017) 1–19

Fig. 14. The z-displacement (deflection) contour for the sandwich plate using proposed method.

Fig. 15b. Convergence study of shear traction T1 of the 2nd interface for m0 , n0 .

Fig. 15a. Convergence study of shear traction T1 of the 1st interface for m0 , n0 .

on the proposed method, the max normal traction stresses for the 1st and 2nd interfaces are 53.297 MPa and 38.418 MPa, respectively. The results from ANSYS and the present method match well.(Fig. 18(a) and (b) Fig. 18(c) and (d))

as: ( ) 𝜉 𝑓𝑇 1 = 𝜉1 1 − 𝜉1 3 ( ) 𝜉 𝑓𝑇 2 = 𝜉2 1 − 𝜉2

6.2.2. Clamp-Free sandwich plate 1

3

𝜉

𝑓𝑇 𝑗 = 1, (𝑖 = 1, 2)(𝑗 = 1, 2) 𝑖

(6–2–4a,b,c) In this Section, the boundary condition functions are defined as:

Thus, the normal traction stress field are compared in Figure 18. It is noted that the max normal traction stresses obtained from ANSYS are 48.0 MPa for the 1st interface and 41.5 MPa for the 2nd interface. Based

𝜉

𝜉

𝜉

𝑓𝑢11 = 𝑓𝑢21 = 𝑓𝑢31 = 1 − 𝜉1 ( ) 𝜉 𝜉 𝜉 𝑓𝑢12 = 𝑓𝑢22 = 𝑓𝑢32 = 𝜉2 1 − 𝜉2

Fig. 16. (a)-(b). The shear traction T1 ’s contour of the 1st and 2nd interface from proposed method. 12

ARTICLE IN PRESS

JID: MS T. Li, R.K. Kapania

[m5GeSdc;August 18, 2017;9:38] International Journal of Mechanical Sciences 000 (2017) 1–19

Fig. 16. (c)-(d). The shear traction T1 ’s contour of the 1st and 2nd interface from ANSYS.

Table 4a Convergence study of displacement about the in-plane displacement order. PO unit:mm ∗

(5,5) 36.077

(6,6) 36.188

(7,7) 37.943

(8,8) 38.164

(9,9) 37.963

(10,10) 37.970

(11,11) 38.171

(a,b) refers to the displacement order in the plate x-y direction.

𝜉

𝜉

𝜉

𝑓𝑢13 = 𝑓𝑢23 = 𝑓𝑢33 = 1

(6–2–5 a,b,c)

Thus, the orthogonal basis function space for the sandwich plate structure will be as follows: 𝜉 𝜑𝑢1𝑖 (𝑖 = 1, 2, 3) is the orthogonal polynomial series with weight function (1 − 𝜉 1 )2 , 𝜉 𝜓𝑢𝑖2 (𝑖 = 1, 2, 3) is the orthogonal polynomial series with weight function 𝜉22 (1 − 𝜉2 )2 , 𝜉

𝜙𝑢1𝑖 (𝑖 = 1, 2, 3) is the orthogonal polynomial series with weight function 1. In this case, the Consistent Orthogonal Basis Functions for 𝜉 2 ,𝜉 3 are still same as those plotted in Figs. 9 and 10. For the basis functions for 𝜉 1 , they are given in Fig. 19. The convergence study for displacement about in-plane displacement order m0 ,n0 is given in Table 4a. The ANSYS reference solution is 38.008 mm. The convergence study for traction T1 and T2 about inplane displacement order m0 ,n0 is given in Table 4b.(Fig. 20(a) and (b) Fig. 21(a) and (b)) The shear traction T1 ,T2 contour is given in Fig. 20–21. The shear traction T1 of x = 0 is plotted against the y-axis in Fig. 22(a) and (b). The shear traction T2 of y = 1/2 is plotted against the x-axis in Fig. 22(c) and (d). From the results in Fig 22 and 23 and Table 4, it is observed that a convergent solution for traction is obtained when the in-plane displacement order is greater than 10.(Fig. 24(a) and (b) Fig. 24(c) and (d)) The normal traction is compared in Fig. 24. The boundary condition functions for traction field are defined as:

Fig. 17. Comparison of shear traction T1 between proposed method and ANSYS. Table 4b Convergence study of shear traction about the in-plane displacement order. PO T1 /MPa T2 /MPa

(7,7) 123.82 73.53

(8,8) 106.60 75.85

(9,9) 114.34 76.34

(10,10) 107.53 75.20

(11,11) 107.23 75.45

6.2.3. Clamp-Free sandwich plate 2 In this Section, the boundary condition functions are defined as: 𝜉

𝜉

𝜉

𝜉

𝜉

𝜉

𝜉

𝑓𝑢11 = 𝑓𝑢21 = 𝑓𝑢31 = 1 − 𝜉1 𝑓𝑢𝜉2 = 𝑓𝑢22 = 𝑓𝑢32 = 𝜉2 1

𝜉

𝑓𝑢13 = 𝑓𝑢23 = 𝑓𝑢33 = 1

(6–2–7 a,b,c)

Thus, the orthogonal basis function space for the sandwich plate structure will be: 𝜉 𝜑𝑢1𝑖 (𝑖 = 1, 2, 3) is the orthogonal polynomial series with weight function (1 − 𝜉 1 )2 , 𝜉 𝜓𝑢𝑖2 (𝑖 = 1, 2, 3) is the orthogonal polynomial series with weight function 2 𝜉2 ,

3

= 1, (𝑖 = 1, 2)(𝑗 = 1, 2)

(6,6) 100.64 67.77

traction for the 1st and 2nd interfaces are 99.1 MPa and 82.7 MPa. The results match well.

( ) 𝜉 𝑓𝑇 1 = 1 − 𝜉1 3 ( ) 𝜉 𝑓𝑇 2 = 𝜉2 1 − 𝜉2 𝜉 𝑓𝑇 𝑗 𝑖

(5,5) 122.82 60.48

(6–2–6 a,b,c)

𝜉

𝜙𝑢1𝑖 (𝑖 = 1, 2, 3) is the orthogonal polynomial series with weight function 1.

The max normal traction for the 1st and 2nd interfaces are 105.48 MPa and 73.542 MPa,respectively. In ANSYS, the max normal 13

JID: MS

ARTICLE IN PRESS

T. Li, R.K. Kapania

[m5GeSdc;August 18, 2017;9:38] International Journal of Mechanical Sciences 000 (2017) 1–19

Fig. 18. (a)-(b). The normal stress contours for the 1st and 2nd interfaces from ANSYS.

Fig. 18. (c)-(d). The normal stress contours for the 1st and 2nd interfaces from the proposed method.

Fig. 19. The Consistent Orthogonal Basis Functions for 𝜉 1 .

Table 5 Convergence study of displacement about the in-plane displacement order.

In this case, the Consistent Orthogonal Basis Functions in terms of 𝜉 1 ,𝜉 2 are still the same as those plotted in Fig. 19. For the basis functions in terms of 𝜉 3 , they are given in Fig. 10. The convergence study of displacement for the in-plane displacement order is given in Table 5. The convergence study of shear traction for in-plane displacement order is given in Table 6. From the result in Table 6, it is observed that the convergence solution for displacement is obtained when the in-plane displacement order is greater than 8. However, for the shear traction solution, the conver-

PO unit:mm

(5,5) 183.94

(6,6) 186.68

(7,7) 188.22

(8,8) 189.57

(9,9) 189.90

(10,10) 190.39

(11,11) 190.59

Table 6 Convergence study of shear traction about the in-plane displacement order. PO T1 /MPa

14

(5,5) 178.79

(6,6) 184.66

(7,7) 189.39

(8,8) 191.98

(9,9) 200.23

(10,10) 203.44

(11,11) 204.86

(12,12) 202.36

JID: MS

ARTICLE IN PRESS

T. Li, R.K. Kapania

[m5GeSdc;August 18, 2017;9:38] International Journal of Mechanical Sciences 000 (2017) 1–19

Fig. 20. (a)-(b). The contour of shear traction T1 for the 1st and 2nd interface.

Fig. 21. (a)-(b). The contour of shear traction T2 for the 1st and 2nd interface.

Fig. 22. (a)-(b). Convergence study of T1 of the 1st and 2nd interface for m0 , n0 .

Fig. 22. (c)-(d). Convergence study of T2 of the 1st and 2nd interface for m0 , n0 . 15

JID: MS

ARTICLE IN PRESS

T. Li, R.K. Kapania

[m5GeSdc;August 18, 2017;9:38] International Journal of Mechanical Sciences 000 (2017) 1–19

Fig. 23. The deflection contour for clamp-free sandwich plate by using the proposed method.

Fig. 24. (a)-(b). The normal stress contours for the 1st and 2nd interfaces from ANSYS.

Fig. 24. (c)-(d). The normal stress contours for the 1st and 2nd interfaces from the proposed method.

gence is obtained when the in-plane displacement order is greater than 11, which is a little larger than that of the last two problems in Section 5.1–5.2. This is because a larger deformation occurs for the case of analysis presented in Section 5.3. The shear traction T1 ’s contour is given in Fig. 25. The shear traction T1 of x = 0 is plotted against the y-axis in Fig. 26. In Fig. 27, the deflection contour is presented.(Fig. 25(a) and (b) Fig. 26(a) and (b))

6.3. Concentrated load In this section, the sandwich plate under concentrated load is studied. The proposed method is a discontinuous finite element, in which there is no concept of nodes. Thus, the application of this method to the case with concentrated forces needs a further discussion. The Dirac’s

16

ARTICLE IN PRESS

JID: MS T. Li, R.K. Kapania

[m5GeSdc;August 18, 2017;9:38] International Journal of Mechanical Sciences 000 (2017) 1–19

Fig. 25. (a)-(b). The contour of shear traction T1 for the 1st and 2nd interface.

Fig. 26. (a)-(b). Convergence study of T1 of the 1st and 2nd interface for m0 , n0 .

Fig. 27. The deflection contour for clamp-free sandwich plate 2 by using the proposed method.

Delta function is used. Thus, the pressure function can be defined as: (

) ( ) 𝑝 = 𝛿 𝑥 − 𝑥0 𝛿 𝑦 − 𝑦0 𝑃0

𝐾𝑄 =

(6-3-1)

where 𝛿(x − x0 ),𝛿(y − y0 ) are the Delta functions located at the (x = x0 ),(y = y0 ), respectively. P0 is the value of the concentrated load. Thus, the load tangential stiffness matrix and the external force vector become: 𝑄=

( ∫𝑡+𝑑𝑡 Γ𝜎

1

1

∫−1 ∫−1

( ( ))( ( )) 𝑝𝑁𝑢𝑇 𝑁𝑛 𝑑 𝜉1 𝑑 𝜉2 = 𝑃0 𝑁𝑢 𝑥0 , 𝑦0 𝑁𝑛 𝑥0 , 𝑦0

(6–3–2 b)

To smooth the solution, a kernel function is applied to approximate the Delta function. Thus, ( ) 𝑃 = 𝐾 𝑥0 , 𝑦0 𝑃0

(

2 2 ( ) − 1 (𝑥−𝑥0 ) +(𝑦−𝑦0 ) 𝐾 𝑥0 , 𝑦0 = 𝑎𝑒 2𝑐 2

)𝑇 (𝑡+𝑑𝑡 ) 𝑡+𝑑𝑡 ( ( ))𝑇 (𝑡+𝑑𝑡 ( )) 𝑛 𝑝𝑑 𝑆 = 𝑁𝑢 𝑥0 , 𝑦0 𝑛 𝑥0 , 𝑦0 𝑁𝑢 𝑡+𝑑𝑡 𝑡+𝑑𝑡

)

(6–3–3 a,b)

where a, c are the parameters √ to control the localization of the kernel function. It is noted that 𝑎𝑐 2𝜋 = 1 must be satisfied such that the in-

(6–3–2 a)

17

ARTICLE IN PRESS

JID: MS T. Li, R.K. Kapania

[m5GeSdc;August 18, 2017;9:38] International Journal of Mechanical Sciences 000 (2017) 1–19

boundary condition of the beam is clamped-free. The free end is loaded by the concentrated force. The transverse deflection of the beam is plotted against x coordinate in Fig. 28. The axial displacement for the top/bottom surface of the layer 3 is plotted against the x coordinate in Fig. 29. The axial displacement is plotted against z coordinate in Fig. 30. The shear stress 𝜏 xz is plotted against z coordinate in Fig. 31. 7. Conclusion The shear traction in laminated plates is analyzed by developing a discontinuous finite element method, in which case the Consistent Orthogonal Basis Functions are developed and applied to represent the displacement functions. The basis functions are very similar to the mode shape functions. The nonlinear finite element analysis with a large deformation is presented for three different boundary conditions. The shear traction is solved accurately by using a much fewer DOFs. The results are also verified with ANSYS. A detailed convergence study is conducted and results are presented. It is observed that the minimum displacement order for a convergent solution of traction is a little larger than that for a convergent displacement solution. However, as compared to the traditional FEM, the number of DOF is much smaller, especially for the convergent traction solution. The present method can be generalized to analyze delaminations of laminated plates. Based on the numerical tests, it is concluded that the proposed basis functions give accurate results even when using much fewer DOFs for large deformation analyses. This is reasonable because the basis functions are very similar to vibration mode shapes, which means that the basis functions derived here consider the fundamental characteristics of the structure (geometry and boundary conditions). In the plots of basis functions, it is observed that the basis functions have larger magnitude in the region where the structure is relatively more flexible. This is consistent with the well-known mechanics principle that the response of

Fig. 28. The transverse deflection of the beam against x coordinate.

tegral of the Kernel function over the effective domain is 1. In this case, the Gaussian Kernel function is applied. It is noted that, if the Kernel function is located at the corner point, the factor a of the Kernel function needs to be multiplied by 4, and if located at the edge (not the corner point), by 2. The problem analyzed in [45,46] is re-solved here. The laminated beam is loaded by the end point concentrated force. The thicknesses for each layers are 2 mm/16 mm/2 mm. The Young’s modulus for each layers are 730 GPa /730 MPa/219 GPa. The shear modulus for each layers are 292 GPa/290 MPa /87.6 GPa. The length of the beam is 100 mm. The

Fig. 29. The axial displacement for the top/bottom surface of the layer 3 against x coordinate.

Fig. 30. The axial displacement against z coordinate.

18

JID: MS

ARTICLE IN PRESS

T. Li, R.K. Kapania

[m5GeSdc;August 18, 2017;9:38] International Journal of Mechanical Sciences 000 (2017) 1–19

Fig. 31. The transverse shear stress against thickness coordinate.

a structure is to minimize the potential energy and the relatively flexible region of the structure always shows a larger magnitude for the displacement response. Thus, the proposed basis functions are able to capture the sharp gradient of the unknown variables by using a much fewer number of DOFs.

[24] Sprague MA, Purkayastha A. Legendre spectral finite elements for Reissner–Mindlin composite plates. Finite Elements Anal Des 2015;105:33–43. [25] Peet YT, Fischer PF. Legendre spectral element method with nearly incompressible materials. Eur J Mech-A/Solids 2014;44:91–103. [26] Rekatsinas CS, Nastos CV, Theodosiou TC, Saravanos DA. A time-domain high-order spectral finite element for the simulation of symmetric and anti-symmetric guided waves in laminated composite strips. Wave Motion 2015;53:1–19. [27] Kant T, Swaminathan K. Analytical solutions for the static analysis of laminated composite and sandwich plates based on a higher order refined theory. Composite Struct 2002;56(4):329–44. [28] Kant T, Swaminathan K. Analytical solutions for free vibration of laminated composite and sandwich plates based on a higher-order refined theory. Composite Struct 2001;53(1):73–85. [29] Pandya BN, Kant T. Higher-order shear deformable theories for flexure of sandwich plates—finite element evaluations. Int J Solids Struct 1988;24(12):1267–86. [30] Reddy JN. A simple higher-order theory for laminated composite plates. J Appl Mech 1984;51(4):745–52. [31] Reddy JN. Analysis of functionally graded plates. Int J Numerical Methods Eng 2000;47(1-3):663–84. [32] Zhang LW, Liew KM, Reddy JN. Postbuckling of carbon nanotube reinforced functionally graded plates with edges elastically restrained against translation and rotation under axial compression. Comput Methods Appl Mech Eng 2016;298:1–28. [33] Carrera E. Historical review of zig-zag theories for multilayered plates and shells. Appl Mech Rev 2003;56(3):287–308. [34] Reddy JN. Mechanics of laminated composite plates and shells: theory and analysis. CRC Press; 2004. [35] Kant T, Swaminathan K. Estimation of transverse/ inter-laminar stresses in laminated composites–a selective review and survey of current developments. Composite Struct 2000;49(1):65–75. [36] Jha DK, Kant T, Singh RK. A critical review of recent research on functionally graded plates. Composite Struct 2013;96:833–49. [37] Li T. On the formulation of a 3-D smooth curved pipe finite element with arbitrary variable cross-section. Thin-Walled Struct 2017;117:317–31. [38] Groh RMJ, Tessler A. Computationally efficient beam elements for accurate stresses in sandwich laminates and laminated composites with delaminations. Comput Methods Appl Mech Eng 2017;320:369–95. [39] Wang X, Shi G. A refined laminated plate theory accounting for the third-order shear deformation and interlaminar transverse stress continuity. Appl Math Model 2015;39(18):5659–80. [40] Pagano NJ. Exact solutions for composite laminates in cylindrical bending. In: Mechanics of composite materials. Netherlands: Springer; 1994. p. 72–85. [41] Pagano NJ. Exact solutions for rectangular bidirectional composites and sandwich plates. In: Mechanics of composite materials. Netherlands: Springer; 1994. p. 86–101. [42] Srinivas S, Rao AK. Bending, vibration and buckling of simply supported thick orthotropic rectangular plates and laminates. Int J Solids Struct 1970;6(11):1463–81. [43] Lee KH, Lin WZ, Chow ST. Bidirectional bending of laminated composite plates using an improved zig-zag model. Composite Struct 1994;28(3):283–94. [44] Mantari JL, Oktem AS, Soares CG. A new higher order shear deformation theory for sandwich and composite laminated plates. Composites Part B 2012;43(3):1489–99. [45] Nallim LG, Oller S, Onate E, Flores FG. A hierarchical finite element element for composite laminated beams using a refined zigzag theory. Composite Struct 2017;163:168–84. [46] Onate E, Wijo A, Oller S. Simple and accurate two-noded beam element for composite laminated beams using a refined zigzag theory. Comput Methods Appl Mech Eng. 2012;213-216:362–82.

Reference [1] Bathe KJ. Finite element procedures. Prentice Hall; 1996. [2] Kapania RK, Raciti S. Recent advances in analysis of laminated beams and plates. Part I-shear effects and buckling. AIAA J 1989;27(7):923–35. [3] Kapania RK, Raciti S. Recent advances in analysis of laminated beams and plates, Part II: vibrations and wave propagation. AIAA J 1989;27(7):935–46. [4] Simo JC, Hughes TJR. Computational inelasticity. Springer; 1998. [5] Goyal VK, Kapania RK. A shear-deformable beam element for the analysis of laminated composites. Finite Elem Anal Des 2007;43(6):463–77. [6] Kapania RK, Yang TY. Buckling, postbuckling, and nonlinear vibrations of imperfect plates. AIAA J 1987;25(10):1338–46. [7] Nosier A, Kapania RK, Reddy JN. Free vibration analysis of laminated plates using a layer wise theory. AIAA J 1993;31(12):2335–46. [8] Zhao W, Kapania RK. Buckling analysis of unitized curvilinearly stiffened composite panels. Composite Struct 2016;135:365–82. [9] Bischoff M, Ramm E. Shear deformable shell elements for large strains and rotations. Int J Numerical Methods Eng 1997;40(23):4427–49. [10] Murakami H. Laminated composite plate theory with improved in-plane responses. J Appl Mech 1986;53(3):661–6. [11] Toledano A, Murakami H. A composite plate theory for arbitrary laminate configurations. J Appl Mech 1987;54(1):181–9. [12] Barbero EJ, Reddy JN. Modeling of delamination in composite laminates using a layer-wise plate theory. Int J Solids Struct 1991;28(3):373–88. [13] Cho M, Parmerter RR. An efficient higher-order plate theory for laminated composites. Composite Struct 1992;20(2):113–23. [14] Reddy JN, Khdeir A. Buckling and vibration of laminated composite plates using various plate theories. AIAA J 1989;27(12):1808–17. [15] Marjanović Miroslav, Djordje Vuksanović. Free vibrations of laminated composite shells using the rotation-free plate elements based on Reddy’s layerwise discontinuous displacement model. Composite Struct 2016;156:320–32. [16] Pramod ALN, Natarajan S, Ferreira AJM, Carrera E, Cinefra M. Static and free vibration analysis of cross-ply laminated plates using the Reissner-mixed variational theorem and the cell based smoothed finite element method. Eur J Mech-A/Solids 2017;62:14–21. [17] McElroy M. Use of an enriched shell finite element to simulate delamination migration in a composite laminate. Composite Struct 2017;167:88–95. [18] Antoine GO, Batra RC. Sensitivity analysis of low-velocity impact response of laminated plates. Int J Impact Eng 2015;78:64–80. [19] Paik S, Gupta SS, Batra RC. Localization of buckling modes in plates and laminates. Composite Struct 2015;120:79–89. [20] Lei ZX, Zhang LW, Liew KM. Analysis of laminated CNT reinforced functionally graded plates using the element-free kp-Ritz method. Composites Part B 2016;84:211–21. [21] Fried I, Malkus DS. Finite element mass matrix lumping by numerical integration with no convergence rate loss. Int J Solids Struct 1975;11(4):461–6. [22] Taylor MA, Wingate BA. A generalized diagonal mass matrix spectral element method for non-quadrilateral elements. Appl Numerical Math 2000;33(1):259–65. [23] Sprague MA, Geers TL. Legendre spectral finite elements for structural dynamics analysis. Commun Numerical Methods Eng 2008;24(12):1953–65.

19