Application of consistent geometric decomposition theorem to dynamic finite element of 3D composite beam based on experimental and numerical analyses

Abstract

Analyzing static and dynamic problems including composite structures has been of high significance in research efforts and industrial applications. In this article, equivalent single layer approach is utilized for dynamic finite element procedures of 3D composite beam as the building block of numerous composite structures. In this model, both displacement and strain fields are decomposed into cross-sectional and longitudinal components, called consistent geometric decomposition theorem. Then, the model is discretized using finite element procedures. Two local coordinate systems and a global one are defined to decouple mechanical degrees of freedom. Furthermore, from the viewpoint of consistent geometric decomposition theorem, the transformation and element mass matrices for those systems are introduced here for the first time. The same decomposition idea can be used for developing element stiffness matrix. Finally, comprehensive validations are conducted for the theory against experimental and numerical results in two case studies and for various conditions.

Keywords

composite beam / dynamic finite element / degrees of freedom coupling / experimental validation / numerical validation

Introduction

Over the preceding years, the use of composite structures for many engineering applications has been growing, raising the necessity for developing and improving appropriate theoretical formulations to model their behaviors [¹–³]. Composite plate [⁴], shell [⁵], and beam [⁶] models have been proposed. On the other hand, dynamic and vibration analyses of structures and mechanisms have been the topic of numerous research attempts for various applications [⁷–¹²] and various mathematical foundations have been developed in solving such engineering problems [¹³–¹⁵]. In addition, flexoelectric and piezoelectric materials and their usages in composite structures have attracted much attention [¹⁶–¹⁹].

Composite beams have attracted much attention as common structural components in mechanical, civil, and aerospace engineering sectors. Static and dynamic analyses based on analytical and approximate methods are documented in the literature for composite beams [²⁰,²¹]. However, it should be noted that using numerical procedures such as finite element method (FEM) is inevitable in many problems including irregular shapes and boundary conditions for which the analytical solutions are not available. Accordingly, the FEM has been frequently utilized to formulate composite beams [²²,²³]. Many composite structures can be thought of being constructed with 3D composite beam elements including a substructure and some patches of different materials or two different layers [²⁴].

Both displacement and strain fields could be formulated using equivalent single layer (ESL) [²⁵,²⁶] or layer-wise [²⁷,²⁸] assumptions each having its own advantages and disadvantages. The accuracy and simplicity of the ESL approach make it important especially for composite structures with less number of layers [²⁹]. A primary composite finite element model based on 3D beam kinematics and ESL assumption was previously derived by Fattahi and Mirdamadi [⁹] and validated with the literature results in which both displacement and strain fields were decomposed into cross-sectional and longitudinal components and then discretized using FEM.

In this article, the above-mentioned model is developed in more details. 3D composite beams are formulated based on Hamilton’s principle for dynamic analyses. At first, an elasticity formulation is developed for various elements. The consistent geometric decomposition theorem is introduced for splitting strain and displacement fields. Then, the continuum model is discretized using FEM and the terms are calculated for a special case; i.e., a 3D beam element. Two local coordinate systems and a global one are introduced to decouple mechanical degrees of freedom and the transformation and mass matrices for those systems are introduced here for the first time. Moreover, the model is validated with experimental and numerical results. The natural frequencies, mode shapes, and responses to base displacements of a cantilevered beam with a piezoelectric layer modeled in ABAQUS software are extracted and compared to those obtained by a MATLAB code based on the theory. In another case study, a cantilevered aluminum beam prototype with two steel patches is made for experimental procedures. A base displacement is applied to the prototype and the tip displacements for a range of excitation frequencies are measured. All numerical and experimental results show a good agreement to those obtained by the theory.

General elastic formulation

In this research, the elastic part of dynamic problems concerning with composite structures are addressed. As discussed in Ref. [⁹], the Hamilton’s principle is used for dynamic analysis of composite structures including a substructure bonded to another layer or some other patches and then, both strain and displacement fields are decomposed into cross-sectional and longitudinal terms using the consistent geometric decomposition theorem, as follows:

(1)\mathit{u}=\overline{\mathit{u}}(y,z)\widehat{\mathit{u}}\left(x\right),\hspace{0.17em}\mathit{S}=\overline{\mathit{S}}(y,z)\widehat{\mathit{S}}\left(x\right),

where \widehat{\mathit{u}} and \widehat{\mathit{S}} are generalized quantities related to displacement and strain, respectively. Also, \overline{\mathit{u}}and \overline{\mathit{S}} are matrices relating the displacement and strain fields to their corresponding generalized quantities, respectively. The structure under study is then discretized using FEM. The displacement field \mathit{u} is related to the displacements of nodes {\mathit{u}}_i using mechanical shape functions {\mathit{N}}_{\mathrm{u}} as expressed in Eq. (2):

(2)\mathit{u}(X,t)={\mathit{N}}_{\mathrm{u}}\left(\mathit{X}\right){\mathit{u}}_i\left(t\right).

The relationships between generalized displacement and generalized strain with the nodal displacement vector are considered as:

(3)\widehat{\mathit{u}}={\mathit{B}}_{\mathrm{d}}{\mathit{u}}_i,\hspace{0.17em}\widehat{\mathit{S}}={\mathit{B}}_{\mathrm{u}}{\mathit{u}}_i.

Depending on the element, the matrices {\mathit{B}}_{\mathrm{d}} and {\mathit{B}}_{\mathrm{u}} have relationships with the shape function {\mathit{N}}_{\mathrm{u}}. By replacing the process in the Hamilton’s principle, the following equations are obtained for element matrices:

(4){\mathit{M}}^{\left(\mathrm{e}\right)}={\displaystyle \underset{V^{\left(\mathrm{e}\right)}}{\int }{\mathit{B}}_{\mathrm{d}}^t{\overline{\mathit{u}}}^t\rho \overline{\mathit{u}}{\mathit{B}}_{\mathrm{d}}\mathrm{d}V}={\displaystyle \underset{L^{\left(\mathrm{e}\right)}}{\int }{\mathit{B}}_{\mathrm{d}}^t\left({\displaystyle \underset{A^{\left(\mathrm{e}\right)}}{\int }{\overline{\mathit{u}}}^t\rho \overline{\mathit{u}}\mathrm{d}A}\right)}{\mathit{B}}_{\mathrm{d}}\mathrm{d}x,

(5){\mathit{K}}^{\left(\mathrm{e}\right)}={\displaystyle \underset{V^{\left(\mathrm{e}\right)}}{\int }{\mathit{B}}_{\mathrm{u}}^t{\overline{\mathit{S}}}^t\mathit{c}\overline{\mathit{S}}{\mathit{B}}_{\mathrm{u}}\mathrm{d}V}={\displaystyle \underset{L^{\left(\mathrm{e}\right)}}{\int }{\mathit{B}}_{\mathrm{u}}^t({\displaystyle \underset{A^{\left(\mathrm{e}\right)}}{\int }{\overline{\mathit{S}}}^t\mathit{c}\overline{\mathit{S}}\mathrm{d}A)}{\mathit{B}}_{\mathrm{u}}}\mathrm{d}x,

(6)f^{\left(\mathrm{e}\right)}={\displaystyle \sum _{i=\mathrm{1}}^{nf}{\mathit{B}}_{\mathrm{d}}^t{\overline{\mathit{u}}}^t\mathit{f}}+{\displaystyle \underset{{\Omega }_{\mathrm{1}}}{\int }{\mathit{B}}_{\mathrm{d}}^t{\overline{\mathit{u}}}^t\mathit{t}\mathrm{d}\Omega }+{\displaystyle \underset{V^{\left(\mathrm{e}\right)}}{\int }{\mathit{B}}_{\mathrm{d}}^t{\overline{\mathit{u}}}^t}\mathit{b}\mathrm{d}V,

where \mathrm{c} is the stiffness tensor, \rho the density, \mathit{M} the mass matrix, \mathit{K} the stiffness matrix, and {\mathit{f}}^{\left(\mathrm{e}\right)} the external mechanical force vector for a set of nfconcentrated forces facting on points with coordinates {\mathit{X}}_i at time t, a surface traction ton the surface \Omega, and a body force b. The superscript (e) stands for element, L the length, A the area of the cross section, and V the volume of the element. Also, the superscript t stands for transpose of a matrix. Finally, the subsequent equation of motion can be obtained by assuming that a base acceleration is applied to the system [⁹]:

(7){{\displaystyle \mathit{M}\ddot{\mathit{u}}}}_i+{{\displaystyle \mathit{C}\dot{\mathit{u}}}}_i+{{\displaystyle \mathit{K}\mathit{u}}}_i=-{\mathit{M}\mathit{T}}_{\mathrm{g}}^t{{\displaystyle \ddot{\mathit{u}}}}_{\mathrm{g}},

\mathit{M} and \mathit{K}are considered to be the assembled mass and stiffness matrices, respectively. {\mathit{T}}_g^tis the transformation matrix [³⁰] and {{\displaystyle \ddot{\mathit{u}}}}_i the vector of acceleration of a reference point on the base. The modal damping theory is utilized and \mathit{C}is the physical damping matrix which is considered to be proportional to mass and stiffness matrices and therefore can become orthogonal using the matrix of mode shapes {\mathit{\Phi }}_{\mathrm{N}}(i.e., \mathit{C}={\mathit{\Phi }}_{\mathrm{N}}^{-t}{\mathit{C}}_{\mathrm{N}}{\mathit{\Phi }}_{\mathrm{N}}^{-\mathrm{1}}); where {\mathit{C}}_{\mathrm{N}}is the diagonal matrix of modal damping. The element on the ith row and ith column of {\mathbf{C}}_{\mathrm{N}}is called C_{\mathrm{N}i}=\mathrm{2}{\zeta }_i{\omega }_iwhich is the modal damping of ith mode; where {\zeta }_iand {\omega }_iare the damping ratio and the natural frequency of ith mode, respectively [³⁰].

Formulation for a 3D composite beam

The decomposition and discretization explained is Section 2 is now performed for a 3D composite beam.

Consistent geometric decomposition theorem

In this paper, two local coordinate systems are used for the cross section of the 3D composite beam shown in Fig. 1:

1) x^{\prime }{y}^{\prime }{z}^{\prime }(the primary local coordinate system) which is the principal one. x^{\prime }is assumed to be the neutral axis of the beam, while y^{\prime } and z^{\prime } are the principal axes of area of the cross section. According to the convention used here, any quantity in this system is distinguished by a ('). Q and C are the neutral point and the shear center of each cross section, respectively, and the axis passing through the shear center of each cross section is denoted by {\overline{x}}^{\prime }.

2) x^{″}{y}^{″}{z}^{″}(the secondary local coordinate system) passing through the center of mass on each cross section (G); x^{″}is parallel to x^{\prime }, and y^{″} and z^{″} are the mass principal axes of inertia of the cross section. According to Fig. 1, the cross section has both geometrical and material symmetries around z^{\prime } axis; accordingly, y^{″}and z^{″} are parallel to y^{\prime } and z^{\prime }, respectively. Based on the convention used, any quantity in the secondary local system is distinguished by a (″).

The Euler-Bernoulli beam and Saint-Venant theories are utilized for bending and torsional effects, respectively. To solve the problem of coupling between axial, bending, and torsional effects, the longitudinal and lateral displacements are considered at O ({{\mathit{u}}^{\prime }}_{\mathrm{o}}) and C ({v^{\prime }}_{\mathrm{C}},{w^{\prime }}_{\mathrm{C}}), respectively; the rotation is also assumed about C ({\theta }_{{\overline{x}}^{\prime }}), and the bending about y^{\prime } and z^{\prime } axes ({\theta }_{y^{\prime }},n{\theta }_{z^{\prime }}). In the primary local coordinate system, the displacement field for the 3D beam is expressed as:

(8){\mathit{u}}^{\prime }=\left\{\begin{array}{l}{\mathit{u}}^{\prime }\hfill \\ v^{\prime }\hfill \\ w^{\prime }\hfill \end{array}\right\}=\left\{\begin{array}{l}{{\mathit{u}}^{\prime }}_{\mathrm{o}}-z^{\prime }{\displaystyle \frac{\partial {w^{\prime }}_{\mathrm{c}}}{\partial x^{\prime }}}+y^{\prime }{\displaystyle \frac{\partial {v^{\prime }}_{\mathrm{c}}}{\partial x^{\prime }}}+\psi {\displaystyle \frac{\partial {\theta }_{{\overline{x}}^{\prime }}}{\partial x^{\prime }}}\hfill \\ \hfill {v^{\prime }}_{\mathrm{c}}-(z^{\prime }-{z^{\prime }}_{\mathrm{c}}){\theta }_{{\overline{x}}^{\prime }}\hfill \\ {w^{\prime }}_{\mathrm{c}}+(y^{\prime }-{y^{\prime }}_{\mathrm{c}}){\theta }_{{\overline{x}}^{\prime }}\hfill \end{array}\right\},

where \psiis the warping function, and {\mathit{u}}^{\prime }, v^{\prime }, and w^{\prime } are displacements in x^{\prime }, y^{\prime }, and z^{\prime } directions, respectively. It should be noted that, all axial, bending and torsional effects considered above must then be transformed to their corresponding parameters in the principal local system using the method explained in Section 3.2. {\sigma }_{x^{\prime }}, {\tau }_{x^{\prime }{z}^{\prime }}, and {\tau }_{x^{\prime }{y}^{\prime }} are the stress components and {\epsilon }_{x^{\prime }}, {\gamma }_{x^{\prime }{z}^{\prime }}, and {\gamma }_{x^{\prime }{y}^{\prime }} are the strain components in the cross section of a 3D composite beam in the principal local coordinate system. Here, the constitutive equations governing on the problem is written as:

(9)\left\{\begin{array}{c}\hfill {\sigma }_{x^{\prime }}\hfill \\ \hfill {\tau }_{x^{\prime }{z}^{\prime }}\hfill \\ \hfill {\tau }_{x^{\prime }{y}^{\prime }}\hfill \end{array}\right\}=\left[\begin{array}{ccc}\hfill E\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill {\mathrm{G}}_{y^{\prime }}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill {\mathrm{G}}_{z^{\prime }}\hfill \end{array}\right]\left\{\begin{array}{c}\hfill {\epsilon }_{x^{\prime }}\hfill \\ \hfill {\gamma }_{x^{\prime }{z}^{\prime }}\hfill \\ \hfill {\gamma }_{x^{\prime }{y}^{\prime }}\hfill \end{array}\right\},

where E is the elasticity modulus and \mathrm{G} the shear modulus. Using small deformation theory and Eq. (8), the vector of strain components \mathit{S}\text{'} can be obtained and decomposed as Ref. [⁹]:

(10){\mathit{S}}^{\prime }=\left\{\begin{array}{c}\hfill {\epsilon }_{x^{\prime }}\hfill \\ \hfill {\gamma }_{x^{\prime }{z}^{\prime }}\hfill \\ \hfill {\gamma }_{x^{\prime }{y}^{\prime }}\hfill \end{array}\right\}=\left\{\begin{array}{c}\hfill {\displaystyle \frac{\partial u^{\prime }}{\partial x^{\prime }}}\hfill \\ \hfill {\displaystyle \frac{\partial {\mathrm{u}}^{\prime }}{\partial z^{\prime }}}+{\displaystyle \frac{\partial w^{\prime }}{\partial x^{\prime }}}\hfill \\ \hfill {\displaystyle \frac{\partial {\mathrm{u}}^{\prime }}{\partial y^{\prime }}}+{\displaystyle \frac{\partial v^{\prime }}{\partial x^{\prime }}}\hfill \end{array}\right\}=\left\{\begin{array}{c}\hfill {\displaystyle \frac{\partial {u^{\prime }}_{\mathrm{o}}}{\partial x^{\prime }}}-z^{\prime }{\displaystyle \frac{{\partial }^{\mathrm{2}}{w^{\prime }}_{\mathrm{c}}}{\partial {x^{\prime }}^{\mathrm{2}}}}+y^{\prime }{\displaystyle \frac{{\partial }^{\mathrm{2}}{v^{\prime }}_{\mathrm{c}}}{\partial {x^{\prime }}^{\mathrm{2}}}}\hfill \\ \hfill \left({\displaystyle \frac{\partial \psi }{\partial z^{\prime }}}+(y^{\prime }-{y^{\prime }}_{\mathrm{c}})\right){\displaystyle \frac{\partial {\theta }_{{\overline{x}}^{\prime }}}{\partial x^{\prime }}}\hfill \\ \hfill \left({\displaystyle \frac{\partial \psi }{\partial y^{\prime }}}+(z^{\prime }-{z^{\prime }}_{\mathrm{c}})\right){\displaystyle \frac{\partial {\theta }_{{\overline{x}}^{\prime }}}{\partial x^{\prime }}}\hfill \end{array}\right\},

(11){\widehat{\mathit{S}}}^{\prime }={\left\{{\displaystyle \frac{\partial {{\mathit{u}}^{\prime }}_{\mathrm{o}}}{\partial x^{\prime }}},-{\displaystyle \frac{{\partial }^{\mathrm{2}}{w^{\prime }}_{\mathrm{c}}}{\partial {x^{\prime }}^{\mathrm{2}}}},{\displaystyle \frac{{\partial }^{\mathrm{2}}{v^{\prime }}_{\mathrm{c}}}{\partial {x^{\prime }}^{\mathrm{2}}}},{\displaystyle \frac{\partial {\theta }_{{\overline{x}}^{\prime }}}{\partial x^{\prime }}}\right\}}^t,

(12){\overline{\mathit{S}}}^{\prime }=\left[\begin{array}{cccc}\hfill \mathrm{1}\hfill & \hfill z^{\prime }\hfill & \hfill -y^{\prime }\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill {\displaystyle \frac{\partial \psi }{\partial z^{\prime }}}+(y^{\prime }-{y^{\prime }}_{\mathrm{c}})\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill {\displaystyle \frac{\partial \psi }{\partial y^{\prime }}}+(z^{\prime }-{z^{\prime }}_{\mathrm{c}})\hfill \end{array}\right].

Similarly, in order to obtain the displacement field in the secondary local system such that the couplings between different degrees of freedom is eliminated, the longitudinal and lateral displacements are considered at G ({{\mathit{u}}^{″}}_{\mathrm{G}}) and C ({v^{″}}_{\mathrm{C}},{w^{″}}_{\mathrm{C}}), respectively. Moreover, the rotation of the cross section is assumed about C ({\theta }_{{\overline{x}}^{″}}), and the bending effects about y^{″} and z^{″} axes passing through G ({\theta }_{y^{″}},{\theta }_{z^{″}}). According to Eq. (8) and by neglecting the effects of warping, the displacement field in the secondary local system is expressed and decomposed as Ref. [⁹]:

(13){\mathit{u}}^{″}=\left\{\begin{array}{c}\hfill {\mathit{u}}^{″}\hfill \\ \hfill v^{″}\hfill \\ \hfill w^{″}\hfill \end{array}\right\}=\left\{\begin{array}{c}\hfill {{\mathit{u}}^{″}}_{\mathrm{G}}-z^{″}{\displaystyle \frac{\partial {w^{″}}_{\mathrm{c}}}{\partial x^{″}}}+y^{″}{\displaystyle \frac{\partial {v^{″}}_{\mathrm{c}}}{\partial x^{″}}}\hfill \\ \hfill {v^{″}}_{\mathrm{c}}-(z^{″}-{z^{″}}_{\mathrm{c}}){\theta }_{{\overline{x}}^{″}}\hfill \\ \hfill {w^{″}}_{\mathrm{c}}+(y^{″}-{y^{″}}_{\mathrm{c}}){\theta }_{{\overline{x}}^{″}}\hfill \end{array}\right\},

(14){\widehat{\mathit{u}}}^{″}\left(x^{″}\right)={\left\{{{\mathit{u}}^{″}}_{\mathrm{G}},{v^{″}}_{\mathrm{c}},{w^{″}}_{\mathrm{c}},{\theta }_{{\overline{x}}^{″}},{\displaystyle \frac{\partial {w^{″}}_{\mathrm{c}}}{\partial x^{″}}},{\displaystyle \frac{\partial {v^{″}}_{\mathrm{c}}}{\partial x^{″}}}\right\}}^t,

(15)\overline{\mathit{u}}(y^{″},z^{″})=\left[\begin{array}{cccccc}\hfill \mathrm{1}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill z^{″}\hfill & \hfill -y^{″}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{1}\hfill & \hfill \mathrm{0}\hfill & \hfill -(z^{″}-{z^{″}}_{\mathrm{c}})\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{1}\hfill & \hfill (y^{″}-{y^{″}}_{\mathrm{c}})\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill \end{array}\right].

The degrees of freedom used in the secondary local coordinate system should also be transformed to their corresponding parameters in the principal local system according to Section 3.2.

Local transformations

From now on, the subsequent formulation is introduced in more details for the first time. To assemble the element matrices, it is vital to state all degrees of freedom in the same local coordinate system. Therefore, the primary degrees of freedom (used in Eq. (8)) and the stiffness matrix must be transformed from shear center to the neutral point:

(16)\{\begin{array}{c}\hfill {v^{\prime }}_{\mathrm{c}}={v^{\prime }}_{\mathrm{o}}-{z^{\prime }}_{\mathrm{c}}{\theta }_{x^{\prime }}，\hfill \\ \hfill {w^{\prime }}_{\mathrm{c}}={w^{\prime }}_{\mathrm{o}}+{y^{\prime }}_{\mathrm{c}}{\theta }_{x^{\prime }}，\hfill \\ \hfill {\theta }_{{\overline{x}}^{\prime }}={\theta }_{x^{\prime }}\mathrm{.}\hfill \end{array}

Accordingly, the following relationship can be written between those degrees of freedom:

{{\mathit{u}}^{\prime }}_{\mathrm{k}}={\mathit{T}}_{\mathrm{k}}{\mathit{u}}^{\prime },

{{\mathit{u}}^{\prime }}_{\mathrm{k}}={\left\{\begin{array}{llllllllllll}{u^{\prime }}_{\mathrm{o}\mathrm{1}}\hfill & {v^{\prime }}_{\mathrm{c}\mathrm{1}}\hfill & {w^{\prime }}_{\mathrm{c}\mathrm{1}}\hfill & {\theta }_{{\overline{x}}^{\prime }\mathrm{1}}\hfill & {\theta }_{y^{\prime }\mathrm{1}}\hfill & {\theta }_{z^{\prime }\mathrm{1}}\hfill & {u^{\prime }}_{\mathrm{o}\mathrm{2}}\hfill & {v^{\prime }}_{\mathrm{c}\mathrm{2}}\hfill & {w^{\prime }}_{\mathrm{c}\mathrm{2}}\hfill & {\theta }_{{\overline{x}}^{\prime }\mathrm{2}}\hfill & {\theta }_{y^{\prime }\mathrm{2}}\hfill & {\theta }_{z^{\prime }\mathrm{2}}\hfill \end{array}\right\}}^t,

{\mathit{u}}^{\prime }={\left\{\begin{array}{llllllllllll}{u^{\prime }}_{\mathrm{o}\mathrm{1}}\hfill & {v^{\prime }}_{\mathrm{o}\mathrm{1}}\hfill & {w^{\prime }}_{\mathrm{o}\mathrm{1}}\hfill & {\theta }_{x^{\prime }\mathrm{1}}\hfill & {\theta }_{y^{\prime }\mathrm{1}}\hfill & {\theta }_{z^{\prime }\mathrm{1}}\hfill & {u^{\prime }}_{\mathrm{o}\mathrm{2}}\hfill & {v^{\prime }}_{\mathrm{o}\mathrm{2}}\hfill & {w^{\prime }}_{\mathrm{o}\mathrm{2}}\hfill & {\theta }_{x^{\prime }\mathrm{2}}\hfill & {\theta }_{y^{\prime }\mathrm{2}}\hfill & {\theta }_{z^{\prime }\mathrm{2}}\hfill \end{array}\right\}}^t,

(17){\mathbf{T}}_{\mathrm{k}}=\left[\begin{array}{cc}\hfill {\mathbf{t}}_{\mathrm{k}}\hfill & \hfill \mathbf{0}\hfill \\ \hfill \mathbf{0}\hfill & \hfill {\mathbf{t}}_{\mathrm{k}}\hfill \end{array}\right]\mathrm{;}\mathrm{}{\mathit{t}}_{\mathrm{k}}=\left[\begin{array}{cccccc}\hfill \mathrm{1}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{1}\hfill & \hfill \mathrm{0}\hfill & \hfill -{z^{\prime }}_{\mathrm{c}}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{1}\hfill & \hfill {y^{\prime }}_{\mathrm{c}}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{1}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{1}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{1}\hfill \end{array}\right].

The stiffness matrix obtained from the primary displacement field ({\mathit{K}}_{\mathrm{e}\mathrm{l}}) is transformed to the stiffness matrix in the principal coordinate system ({{\mathit{K}}^{\prime }}^{\left(\mathrm{e}\right)}):

(18){{\mathit{K}}^{\prime }}^{\left(\mathrm{e}\right)}={\mathit{T}}_{\mathrm{k}}^t{\mathit{K}}_{\mathrm{e}\mathrm{l}}{\mathit{T}}_{\mathrm{k}}.

Similarly, the secondary degrees of freedom (used in Eq. (13)) and the mass matrix can be transformed from the mass and shear centers to the neutral point:

(19)\{\begin{array}{c}\hfill {{\mathit{u}}^{″}}_{\mathrm{G}}={{\mathit{u}}^{\prime }}_{\mathrm{G}}={{\mathit{u}}^{\prime }}_{\mathrm{o}}+{z^{\prime }}_{\mathrm{G}}{\theta }_{y^{\prime }}-{y^{\prime }}_{\mathrm{G}}{\theta }_{z^{\prime }},\hfill \\ \hfill {v^{″}}_{\mathrm{c}}={v^{\prime }}_{\mathrm{c}}={v^{\prime }}_{\mathrm{o}}-{z^{\prime }}_{\mathrm{c}}{\theta }_{x^{\prime }},\hfill \\ \hfill {w^{″}}_{\mathrm{c}}={w^{\prime }}_{\mathrm{c}}={w^{\prime }}_{\mathrm{o}}+{y^{\prime }}_{\mathrm{c}}{\theta }_{x^{\prime }},\hfill \\ \hfill {\theta }_{{\overline{x}}^{\prime }}={\theta }_{x^{\prime }},\hfill \\ \hfill {\theta }_{y^{″}}={\theta }_{y^{\prime }},\hfill \\ \hfill {\theta }_{z^{″}}={\theta }_{z^{\prime }}.\hfill \end{array}

Hence, the relationship between those degrees of freedom can be expressed in the matrix form as:

{{\mathit{u}}^{\prime }}_{\mathrm{m}}={\mathit{T}}_{\mathrm{m}}{\mathit{u}}^{\prime },

{{\mathit{u}}^{\prime }}_{\mathrm{m}}={\left\{\begin{array}{llllllllllll}{u^{″}}_{\mathrm{G}\mathrm{1}}\hfill & {v^{″}}_{\mathrm{c}\mathrm{1}}\hfill & {w^{″}}_{\mathrm{c}\mathrm{1}}\hfill & {\theta }_{{\overline{x}}^{″}\mathrm{1}}\hfill & {\theta }_{y^{″}\mathrm{1}}\hfill & {\theta }_{z^{″}\mathrm{1}}\hfill & {{\mathit{u}}^{″}}_{\mathrm{G}\mathrm{2}}\hfill & {v^{″}}_{\mathrm{c}\mathrm{2}}\hfill & {w^{″}}_{\mathrm{c}\mathrm{2}}\hfill & {\theta }_{{\overline{x}}^{″}\mathrm{2}}\hfill & {\theta }_{y^{″}\mathrm{2}}\hfill & {\theta }_{z^{″}\mathrm{2}}\hfill \end{array}\right\}}^t,

{\mathit{T}}_{\mathrm{m}}=\left[\begin{array}{cc}\hfill {\mathit{t}}_{\mathrm{m}}\hfill & \hfill \mathit{0}\hfill \\ \hfill \mathit{0}\hfill & \hfill {\mathit{t}}_{\mathrm{m}}\hfill \end{array}\right]\mathrm{;}

(20){\mathit{t}}_{\mathrm{m}}=\left[\begin{array}{cccccc}\hfill \mathrm{1}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill {z^{\prime }}_{\mathrm{G}}\hfill & \hfill -{y^{\prime }}_{\mathrm{G}}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{1}\hfill & \hfill \mathrm{0}\hfill & \hfill -{z^{\prime }}_{\mathrm{c}}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{1}\hfill & \hfill {y^{\prime }}_{\mathrm{c}}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{1}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{1}\hfill & \hfill \mathrm{0}\hfill \\ \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{0}\hfill & \hfill \mathrm{1}\hfill \end{array}\right].

Lumping procedure

To discretize, composite 3D beam elements with two nodes and six degrees of freedom per each node are used. The vector of mechanical degrees of freedom in the principal local coordinate system is written as Refs. [³⁰,³¹]:

(21)\underset{\mathrm{12}\times \mathrm{1}}{{{\mathit{u}}^{\prime }}_i}\equiv {\left[\begin{array}{llllllllllll}{{\mathit{u}}^{\prime }}_{\mathrm{1}}\hfill & {v^{\prime }}_{\mathrm{1}}\hfill & {w^{\prime }}_{\mathrm{1}}\hfill & {{\theta }^{\prime }}_{x\mathrm{1}}\hfill & {{\theta }^{\prime }}_{y\mathrm{1}}\hfill & {{\theta }^{\prime }}_{z\mathrm{1}}\hfill & {u^{\prime }}_{\mathrm{2}}\hfill & {v^{\prime }}_{\mathrm{2}}\hfill & {w^{\prime }}_{\mathrm{2}}\hfill & {{\theta }^{\prime }}_{x\mathrm{2}}\hfill & {{\theta }^{\prime }}_{y\mathrm{2}}\hfill & {{\theta }^{\prime }}_{z\mathrm{2}}\hfill \end{array}\right]}^t,

where the parameters \mathrm{u}, v, and ware displacements in the x, y, and z directions, respectively. Furthermore, {\theta }_x, {\theta }_y, and {\theta }_z are rotations about x, y, and z directions, respectively; and the subscripts 1 and 2 represent the number of node. Six functions (three displacements and three rotations) are generally essential for analyzing the kinematics of such an element; however, four principal functions (u, v, w, and {\theta }_x) are sufficient here because the Euler-Bernoulli beam theory is utilized. Linear shape functions are usually appropriate for u and {\theta }_x, while cubic shape functions are mostly adopted for lateral displacements v and w [²³,²⁴]. Now, the displacement field is discretized as:

(22)\underset{\mathrm{4}\times \mathrm{1}}{{\mathit{u}}^{\prime }}(x^{\prime },t)\equiv \left\{\begin{array}{l}{u}^{\prime }(x^{\prime },t)\hfill \\ v^{\prime }(x^{\prime },t)\hfill \\ w^{\prime }(x^{\prime },t)\hfill \\ {{\theta }^{\prime }}_x(x^{\prime },t)\hfill \end{array}\right\}={{\mathit{N}}^{\prime }}_{\mathrm{u}}\left(x^{\prime }\right){{\mathit{u}}^{\prime }}_i\left(t\right).

The linear and cubic shape functions are presented in Eqs. (23) and (24), respectively:

(23)N_1^{\prime 1}=1-\frac{x^{\prime }}{L},N_2^{\prime 1}=1-\frac{x^{\prime }}{L^{\prime }}

(24){N^{\prime }}_{\mathrm{3}}^c={\displaystyle \frac{\mathrm{1}}{L^{\mathrm{3}}}}(-\mathrm{2}{x^{\prime }}^{\mathrm{3}}+\mathrm{3}{x^{\prime }}^{\mathrm{2}}L),{N^{\prime }}_{\mathrm{4}}^c={\displaystyle \frac{\mathrm{1}}{L^{\mathrm{3}}}}({x^{\prime }}^{\mathrm{3}}L-{x^{\prime }}^{\mathrm{2}}{L}^{\mathrm{2}}),

The matrix of shape functions {{\mathit{N}}^{\prime }}_{\mathrm{u}} in Eq. (25) is expressed in the following form:

(25)\underset{\mathrm{4}\times \mathrm{12}}{{N^{\prime }}_u}=\left[\begin{array}{llllllllllll}{N^{\prime }}_{\mathrm{1}}^l& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& {N^{\prime }}_{\mathrm{2}}^l& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {N^{\prime }}_{\mathrm{1}}^c& \mathrm{0}& \mathrm{0}& \mathrm{0}& {N^{\prime }}_{\mathrm{2}}^c& \mathrm{0}& {N^{\prime }}_{\mathrm{3}}^c& \mathrm{0}& \mathrm{0}& \mathrm{0}& {N^{\prime }}_{\mathrm{4}}^c\\ \mathrm{0}& \mathrm{0}& {N^{\prime }}_{\mathrm{1}}^c& \mathrm{0}& -{N^{\prime }}_{\mathrm{2}}^c& \mathrm{0}& \mathrm{0}& \mathrm{0}& {N^{\prime }}_{\mathrm{3}}^c& \mathrm{0}& -{N^{\prime }}_{\mathrm{4}}^c& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& {N^{\prime }}_{\mathrm{1}}^l& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& {N^{\prime }}_{\mathrm{2}}^l& \mathrm{0}& \mathrm{0}\end{array}\right],

{{\mathit{B}}^{\prime }}_{\mathrm{u}} and {{\mathit{B}}^{\prime }}_{\mathrm{d}} can be obtained using Eqs. (3, 11, 21) and (3, 14, 21), respectively:

(26)\underset{\mathrm{4}\times \mathrm{12}}{{B^{\prime }}_u}=\left[\begin{array}{llllllllllll}{N^{\prime }}_{\mathrm{1},x^{\prime }}^l& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& {N^{\prime }}_{\mathrm{2},x^{\prime }}^l& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -{N^{\prime }}_{\mathrm{1},x^{\prime }{x}^{\prime }}^c& \mathrm{0}& {N^{\prime }}_{\mathrm{2},x^{\prime }{x}^{\prime }}^c& \mathrm{0}& \mathrm{0}& \mathrm{0}& -{N^{\prime }}_{\mathrm{3},x^{\prime }{x}^{\prime }}^c& \mathrm{0}& {N^{\prime }}_{\mathrm{4},x^{\prime }{x}^{\prime }}^c& \mathrm{0}\\ \mathrm{0}& {N^{\prime }}_{\mathrm{1},x^{\prime }{x}^{\prime }}^c& \mathrm{0}& \mathrm{0}& \mathrm{0}& {N^{\prime }}_{\mathrm{2},x^{\prime }{x}^{\prime }}^c& \mathrm{0}& {N^{\prime }}_{\mathrm{3},x^{\prime }{x}^{\prime }}^c& \mathrm{0}& \mathrm{0}& \mathrm{0}& {N^{\prime }}_{\mathrm{4},x^{\prime }{x}^{\prime }}^c\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& {N^{\prime }}_{\mathrm{1},x^{\prime }}^l& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& {N^{\prime }}_{\mathrm{2},x^{\prime }}^l& \mathrm{0}& \mathrm{0}\end{array}\right],

(27){B^{\prime }}_{\mathrm{d}}=\left[\begin{array}{llllllllllll}{N^{\prime }}_{\mathrm{1}}^l& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& {N^{\prime }}_{\mathrm{2}}^l& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {N^{\prime }}_{\mathrm{1}}^c& \mathrm{0}& \mathrm{0}& \mathrm{0}& {N^{\prime }}_{\mathrm{2}}^c& \mathrm{0}& {N^{\prime }}_{\mathrm{3}}^c& \mathrm{0}& \mathrm{0}& \mathrm{0}& {N^{\prime }}_{\mathrm{4}}^c\\ \mathrm{0}& \mathrm{0}& {N^{\prime }}_{\mathrm{1}}^c& \mathrm{0}& -{N^{\prime }}_{\mathrm{2}}^c& \mathrm{0}& \mathrm{0}& \mathrm{0}& {N^{\prime }}_{\mathrm{3}}^c& \mathrm{0}& -{N^{\prime }}_{\mathrm{4}}^c& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& {N^{\prime }}_{\mathrm{1}}^l& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& {N^{\prime }}_{\mathrm{2}}^l& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& -{N^{\prime }}_{\mathrm{1},x^{\prime }}^c& \mathrm{0}& {N^{\prime }}_{\mathrm{2},x^{\prime }}^c& \mathrm{0}& \mathrm{0}& \mathrm{0}& -{N^{\prime }}_{\mathrm{3},x^{\prime }}^c& \mathrm{0}& {N^{\prime }}_{\mathrm{4},x^{\prime }}^c& \mathrm{0}\\ \mathrm{0}& {N^{\prime }}_{\mathrm{1},x^{\prime }}^c& \mathrm{0}& \mathrm{0}& \mathrm{0}& {N^{\prime }}_{\mathrm{2},x^{\prime }}^c& \mathrm{0}& {N^{\prime }}_{\mathrm{3},x^{\prime }}^c& \mathrm{0}& \mathrm{0}& \mathrm{0}& {N^{\prime }}_{\mathrm{4},x^{\prime }}^c\end{array}\right].

All parameters needed to calculate the element stiffness and mass matrices are now available. For the primary displacement field, the element mass matrix is obtained in Eqs. (28) and (29) using Eq. (4), which should then be transformed to the principal local and then the global coordinate systems as explained in Section 3.2. It should be noted that the stiffness matrix can be similarly obtained.

(28){\mathit{M}}_{\mathrm{e}\mathrm{l}}={\displaystyle \frac{\mathrm{1}}{L^{{\left(\mathrm{e}\right)}^{\mathrm{3}}}}}\left[\begin{array}{ll}{{\mathit{M}}^{\prime }}_{\mathrm{11}}\hfill & {{\mathit{M}}^{\prime }}_{\mathrm{12}}\hfill \\ {{\mathit{M}}^{\prime }}_{\mathrm{21}}\hfill & {{\mathit{M}}^{\prime }}_{\mathrm{22}}\hfill \end{array}\right],

{M^{\prime }}_{\mathrm{11}}=\left[\begin{array}{llllll}{\displaystyle \frac{{{\rm P}}_m{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{3}}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {\displaystyle \frac{\mathrm{13}{{\rm P}}_m{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{35}}}+{\displaystyle \frac{\mathrm{6}{{\rm P}}_{z^{″}}{L}^{{(e)}^{\mathrm{2}}}}{\mathrm{5}}}& \mathrm{0}& {\displaystyle \frac{\mathrm{7}{{\rm P}}_m{z^{″}}_c{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{20}}}& \mathrm{0}& {\displaystyle \frac{\mathrm{11}{{\rm P}}_m{L}^{{(e)}^{\mathrm{5}}}}{\mathrm{210}}}+{\displaystyle \frac{{{\rm P}}_{z^{″}}{L}^{{(e)}^{\mathrm{3}}}}{\mathrm{10}}}\\ \mathrm{0}& \mathrm{0}& {\displaystyle \frac{\mathrm{13}{{\rm P}}_m{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{70}}}+{\displaystyle \frac{\mathrm{6}{{\rm P}}_{y^{″}}{L}^{{(e)}^{\mathrm{2}}}}{\mathrm{5}}}& -{\displaystyle \frac{\mathrm{7}{{\rm P}}_m{y^{″}}_c{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{20}}}& -{\displaystyle \frac{\mathrm{11}{{\rm P}}_m{L}^{{(e)}^{\mathrm{5}}}}{\mathrm{210}}}-{\displaystyle \frac{{{\rm P}}_{y^{″}}{L}^{{(e)}^{\mathrm{3}}}}{\mathrm{10}}}& \mathrm{0}\\ \mathrm{0}& {\displaystyle \frac{\mathrm{7}{{\rm P}}_m{z^{″}}_c{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{20}}}& -{\displaystyle \frac{\mathrm{7}{{\rm P}}_m{y^{″}}_c{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{20}}}& {\displaystyle \frac{{{\rm P}}_t{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{3}}}& {\displaystyle \frac{{{\rm P}}_m{y^{″}}_c{L}^{{(e)}^{\mathrm{5}}}}{\mathrm{20}}}& {\displaystyle \frac{{{\rm P}}_m{z^{″}}_c{L}^{{(e)}^{\mathrm{5}}}}{\mathrm{20}}}\\ \mathrm{0}& \mathrm{0}& -{\displaystyle \frac{\mathrm{11}{{\rm P}}_m{L}^{{(e)}^{\mathrm{5}}}}{\mathrm{210}}}-{\displaystyle \frac{{{\rm P}}_{y^{″}}{L}^{{(e)}^{\mathrm{3}}}}{\mathrm{10}}}& {\displaystyle \frac{{{\rm P}}_m{y^{″}}_c{L}^{{(e)}^{\mathrm{5}}}}{\mathrm{20}}}& {\displaystyle \frac{{{\rm P}}_m{L}^{{(e)}^{\mathrm{6}}}}{\mathrm{105}}}+{\displaystyle \frac{\mathrm{2}{{\rm P}}_{y^{″}}{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{15}}}& \mathrm{0}\\ \mathrm{0}& {\displaystyle \frac{\mathrm{11}{{\rm P}}_m{L}^{{(e)}^{\mathrm{5}}}}{\mathrm{210}}}+{\displaystyle \frac{{{\rm P}}_{z^{″}}{L}^{{(e)}^{\mathrm{3}}}}{\mathrm{10}}}& \mathrm{0}& {\displaystyle \frac{{{\rm P}}_m{z^{″}}_c{L}^{{(e)}^{\mathrm{5}}}}{\mathrm{20}}}& \mathrm{0}& {\displaystyle \frac{{{\rm P}}_m{L}^{{(e)}^{\mathrm{6}}}}{\mathrm{105}}}+{\displaystyle \frac{\mathrm{2}{{\rm P}}_{z^{″}}{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{15}}}\end{array}\right],

{M^{\prime }}_{\mathrm{12}}=\left[\begin{array}{llllll}{\displaystyle \frac{{{\rm P}}_m{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{6}}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {\displaystyle \frac{\mathrm{9}{{\rm P}}_m{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{70}}}-{\displaystyle \frac{\mathrm{6}{{\rm P}}_{z^{″}}{L}^{{(e)}^{\mathrm{2}}}}{\mathrm{5}}}& \mathrm{0}& {\displaystyle \frac{\mathrm{3}{{\rm P}}_m{z^{″}}_c{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{20}}}& \mathrm{0}& {\displaystyle \frac{\mathrm{13}{{\rm P}}_m{L}^{{(e)}^{\mathrm{5}}}}{\mathrm{420}}}-{\displaystyle \frac{{{\rm P}}_{z^{″}}{L}^{{(e)}^{\mathrm{3}}}}{\mathrm{10}}}\\ \mathrm{0}& \mathrm{0}& {\displaystyle \frac{\mathrm{9}{{\rm P}}_m{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{70}}}-{\displaystyle \frac{\mathrm{6}{{\rm P}}_{y^{″}}{L}^{{(e)}^{\mathrm{2}}}}{\mathrm{5}}}& -{\displaystyle \frac{\mathrm{3}{{\rm P}}_m{y^{″}}_c{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{20}}}& -{\displaystyle \frac{\mathrm{13}{{\rm P}}_m{L}^{{(e)}^{\mathrm{5}}}}{\mathrm{420}}}+{\displaystyle \frac{{{\rm P}}_{y^{″}}{L}^{{(e)}^{\mathrm{3}}}}{\mathrm{10}}}& \mathrm{0}\\ \mathrm{0}& {\displaystyle \frac{\mathrm{3}{{\rm P}}_m{z^{″}}_c{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{20}}}& -{\displaystyle \frac{\mathrm{3}{{\rm P}}_m{y^{″}}_c{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{20}}}& {\displaystyle \frac{{{\rm P}}_t{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{6}}}& {\displaystyle \frac{{{\rm P}}_m{y^{″}}_c{L}^{{(e)}^{\mathrm{5}}}}{\mathrm{30}}}& {\displaystyle \frac{{{\rm P}}_m{z^{″}}_c{L}^{{(e)}^{\mathrm{5}}}}{\mathrm{30}}}\\ \mathrm{0}& \mathrm{0}& {\displaystyle \frac{\mathrm{13}{{\rm P}}_m{L}^{{(e)}^{\mathrm{5}}}}{\mathrm{420}}}-{\displaystyle \frac{{{\rm P}}_{y^{″}}{L}^{{(e)}^{\mathrm{3}}}}{\mathrm{10}}}& -{\displaystyle \frac{{{\rm P}}_m{y^{″}}_c{L}^{{(e)}^{\mathrm{5}}}}{\mathrm{30}}}& -{\displaystyle \frac{{{\rm P}}_m{L}^{{(e)}^{\mathrm{6}}}}{\mathrm{140}}}+{\displaystyle \frac{{{\rm P}}_{y^{″}}{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{30}}}& \mathrm{0}\\ \mathrm{0}& {\displaystyle \frac{\mathrm{13}{{\rm P}}_m{L}^{{(e)}^{\mathrm{5}}}}{\mathrm{420}}}+{\displaystyle \frac{{{\rm P}}_{z^{″}}{L}^{{(e)}^{\mathrm{3}}}}{\mathrm{10}}}& \mathrm{0}& -{\displaystyle \frac{{{\rm P}}_m{z^{″}}_c{L}^{{(e)}^{\mathrm{5}}}}{\mathrm{30}}}& \mathrm{0}& -{\displaystyle \frac{{{\rm P}}_m{L}^{{(e)}^{\mathrm{6}}}}{\mathrm{140}}}+{\displaystyle \frac{{{\rm P}}_{z^{″}}{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{30}}}\end{array}\right],

(29){M^{\prime }}_{\mathrm{22}}=\left[\begin{array}{llllll}{\displaystyle \frac{{{\rm P}}_m{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{3}}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {\displaystyle \frac{\mathrm{13}{{\rm P}}_m{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{35}}}+{\displaystyle \frac{\mathrm{6}{{\rm P}}_{z^{″}}{L}^{{(e)}^{\mathrm{2}}}}{\mathrm{5}}}& \mathrm{0}& {\displaystyle \frac{\mathrm{7}{{\rm P}}_m{z^{″}}_c{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{20}}}& \mathrm{0}& -{\displaystyle \frac{\mathrm{11}{{\rm P}}_m{L}^{{(e)}^{\mathrm{5}}}}{\mathrm{210}}}-{\displaystyle \frac{{{\rm P}}_{z^{″}}{L}^{{(e)}^{\mathrm{3}}}}{\mathrm{10}}}\\ \mathrm{0}& \mathrm{0}& {\displaystyle \frac{\mathrm{13}{{\rm P}}_m{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{35}}}+{\displaystyle \frac{\mathrm{6}{{\rm P}}_{y^{″}}{L}^{{(e)}^{\mathrm{2}}}}{\mathrm{5}}}& -{\displaystyle \frac{\mathrm{7}{{\rm P}}_m{y}^{″}{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{20}}}& {\displaystyle \frac{\mathrm{11}{{\rm P}}_m{L}^{{(e)}^{\mathrm{5}}}}{\mathrm{210}}}+{\displaystyle \frac{{{\rm P}}_{y^{″}}{L}^{{(e)}^{\mathrm{3}}}}{\mathrm{10}}}& \mathrm{0}\\ \mathrm{0}& {\displaystyle \frac{\mathrm{7}{{\rm P}}_m{z^{″}}_c{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{20}}}& -{\displaystyle \frac{\mathrm{7}{{\rm P}}_m{y}^{″}{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{20}}}& {\displaystyle \frac{{{\rm P}}_t{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{3}}}& -{\displaystyle \frac{{{\rm P}}_m{y^{″}}_c{L}^{{(e)}^{\mathrm{5}}}}{\mathrm{20}}}& -{\displaystyle \frac{{{\rm P}}_m{z^{″}}_c{L}^{{(e)}^{\mathrm{5}}}}{\mathrm{20}}}\\ \mathrm{0}& \mathrm{0}& {\displaystyle \frac{\mathrm{11}{{\rm P}}_m{L}^{{(e)}^{\mathrm{5}}}}{\mathrm{210}}}+{\displaystyle \frac{{{\rm P}}_{y^{″}}{L}^{{(e)}^{\mathrm{3}}}}{\mathrm{10}}}& -{\displaystyle \frac{{{\rm P}}_m{y^{″}}_c{L}^{{(e)}^{\mathrm{5}}}}{\mathrm{20}}}& {\displaystyle \frac{{{\rm P}}_m{L}^{{(e)}^{\mathrm{6}}}}{\mathrm{105}}}+{\displaystyle \frac{\mathrm{2}{{\rm P}}_{y^{″}}{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{15}}}& \mathrm{0}\\ \mathrm{0}& -{\displaystyle \frac{\mathrm{11}{{\rm P}}_m{L}^{{(e)}^{\mathrm{5}}}}{\mathrm{210}}}-{\displaystyle \frac{{{\rm P}}_{z^{″}}{L}^{{(e)}^{\mathrm{3}}}}{\mathrm{10}}}& \mathrm{0}& -{\displaystyle \frac{{{\rm P}}_m{z^{″}}_c{L}^{{(e)}^{\mathrm{5}}}}{\mathrm{20}}}& \mathrm{0}& {\displaystyle \frac{{{\rm P}}_m{L}^{{(e)}^{\mathrm{6}}}}{\mathrm{105}}}+{\displaystyle \frac{\mathrm{2}{{\rm P}}_{z^{″}}{L}^{{(e)}^{\mathrm{4}}}}{\mathrm{15}}}\end{array}\right],

where P_{\mathrm{m}}, P_{y^{″}}, P_{z^{″}}, and P_t are obtained based on Eq. (4) using the following equations:

(30)P_{\mathrm{m}}={\displaystyle \underset{A}{\int }\rho \mathrm{d}A},

(31)P_{y^{″}}={\displaystyle \underset{A}{\int }\rho {z^{″}}^{\mathrm{2}}\mathrm{d}A},

(32)P_{z^{″}}={\displaystyle \underset{A}{\int }\rho {y^{″}}^{\mathrm{2}}\mathrm{d}A},

(33)P_t={\displaystyle \underset{A}{\int }\rho ({(z^{″}-{z^{″}}_{\mathrm{c}})}^{\mathrm{2}}+{(y^{″}-{y^{″}}_{\mathrm{c}})}^{\mathrm{2}})\mathrm{d}A}.

Results and discussion

In this section, two case studies are performed in order to validate the proposed theory with ABAQUS and experimental results.

Numerical validations

Composite structures composed of a substructure and some piezoelectric patches or layers are good examples to be analyzed using the proposed theory. In many smart structures including piezoelectric patches or layers, vibration analysis is typically a required and appealing trend. Predicting and extracting natural frequencies and mode shapes of the structure is actually an indispensable and substantial step toward analyzing and designing energy harvesters and for other applications [⁷,²⁴,³²]. As is common in the literature, electrical and electromechanical characteristics of piezoelectric materials can usually be neglected when extracting natural frequencies and mode shapes which is equivalent to a short circuit electrical boundary condition [³³]. It is an acceptable assumption due to the fact that these characteristics have very small effects on these parameters especially when the thickness of piezoelectric material is small compared to that of sub-structure [³⁴,³⁵]. For instance, in Ref. [³⁵], it is discussed with various examples that by considering the electrical and coupling features, the natural frequencies have slight changes even between short circuit and open circuit conditions. The open circuit natural frequency is slightly higher than that of short circuit and the open circuit frequency actually increases the stiffness of the structure. Hence, the theory developed in the present paper can be a powerful means for extracting natural frequencies and mode shapes of piezoelectric smart structures with an acceptable accuracy.

A composite cantilevered beam including a substructure, a piezoelectric layer, and a tip mass is modeled in ABAQUS software (Fig. 2). The specifications of the structure are presented in Table 1.

The first to sixth un-damped natural frequencies obtained from the MATLAB code based on the proposed formulation are compared to those of ABAQUS. The results are demonstrated graphically in Fig. 3(a). “Bending y” and “Bending z” represent bending modes about y and z axes, respectively. The error between theoretical model and ABAQUS is illustrated in Fig. 3(b).

From Fig. 3, it is obvious that the calculated natural frequencies using the developed model overlap with those of ABAQUS results for bending modes. As explained in Section 3.1, the error in torsional mode (mode 6) was expected to be comparatively higher than the other modes because of neglecting the warping function.

The mode shapes 1–3 of the composite beam under investigation are also shown and compared in Figs. 4–6, respectively. According to Figs. 4–6, the mode shapes calculated from theory are in good agreement with those obtained using ABAQUS software. It should be noted that the ABAQUS model is meshed using solid elements.

In the next step, a harmonic base displacement with an amplitude equal to 0.01 mm is applied to the model and the results are compared with an ABAQUS simulation in three conditions. Modal damping theory with a damping ratio of 0.03 is used for both the aluminum substructure [²³,³⁶] and the piezoelectric layer [³⁵]. The same specifications as Table 1 are used for the model for three different thicknesses of the piezoelectric layer. Another point to be mentioned here is that in the case of applying long-term cyclic loading to piezoceramics as brittle materials damages to to fatigue may be of importance and worthy of investigation. This problem has attracted enormous research attempts [³⁷–⁴⁰] which is not required for and not the focus of the present study. Furthermore, many piezo-polymers and piezo-composites have been manufactured with higher flexibility and less brittleness. The results of (a_{z_{\mathrm{rel}}}/a_{z_{\mathrm{base}}}) over a frequency range are depicted in Figs. 7(a)–7(c) for c = 0.1 mm, c = 0.2 mm, and c = 0.4 mm, respectively; in which a_{z_{\mathrm{rel}}}=a_{z_{\mathrm{tip}}}-a_{z_{\mathrm{base}}}, a_{z_{\mathrm{base}}} is the acceleration of a point on the base in z-direction and a_{z_{\mathrm{t}\mathrm{i}\mathrm{p}}} the acceleration of a point on the tip mass in z-direction. As is obvious from Fig. 7, the results of the proposed theory follow those obtained from ABAQUS with small amounts of error.

For more parametric sudies and numerical validations, the same base excitation is applied to the beam for three different active lengths: L = 15 mm, L = 24 mm, and L = 33 mm and the results are demonstrated in Fig. 8. c = 0.1 mm is considered for all cases.

Moreover, for two different materials of the substructure (aluminum and steel), L = 33 mm, c = 0.1 mm and the same excitation, the results of ABAQUS and MATLAB code (the theoretical model) are depicted in Fig. 9. These results also imply that the theoretical results are in good agreement with those of ABAQUS model.

Experimental validations

In another case study, a cantilevered aluminum beam prototype with a steel tip mass is made. Two identical steel patches are completely bonded to the upper surface of the beam in specified locations (Fig. 10). The specifications of the structure are shown in Table 2.

The prototype is then mounted on a vibrator type VP.4 using a clamped fixture. The whole test setup is demonstrated in Fig. 11. A harmonic base displacement with three different amplitudes is applied to the structure in a frequency range of 0–60 Hz. A damping ratio of 0.03 is considered for both the aluminum substructure and the steel patches [³⁰,³⁶]. An MPU-6050 sensor is bonded on the tip mass to measure the acceleration. The results of tip acceleration are then integrated twice to obtain tip displacement.

For excitation amplitudes of 0.07, 0.1, and 0.13 mm, the results of tip displacement obtained by the proposed theory are compared to those measured in the experimental procedure (Fig. 12). Based on the figure, the errors between theory and experiment are acceptable showing that the theory is in a good agreement with experimental results.

Conclusions

The consistent geometric decomposition theorem was applied to a general dynamic finite element formulation. Afterwards, the continuum model was discretized using FEM for a 3D beam element. The fundamental formulation was previously derived based on ESL approach and validated against literature results. In this paper, the emphasis was on detailed developing and experimentally and numerically validating the theory. Furthermore, the transformation and mass matrices were introduced here for the first time. The natural frequencies and mode shapes of bending modes obtained by the ABAQUS model using solid elements were in good agreement with those calculated by the proposed theory. Actually, only few elements can be utilized to model these structures with a simple MATLAB code and high accuracy rather than using a large number of tetrahedral or quadrilateral solid elements available in commercial finite element packages. As was discussed, the error in torsional modes were due to neglecting the warping function which is suggested to be taken into consideration in future works. Then, numerical validations and various case studies on thickness, length, and material of a composite cantilever beam with a piezoelectric layer was performed for a base acceleration and the results were compared to those of ABAQUS showing good agreement. Moreover, in another case study, a composite beam with two steel patches was applied to base excitations and the results of theory were validated against those obtained by experimental procedures.

References

[1] Madeira J F A, Araújo A L, Mota Soares C M, Mota Soares C A. Multiobjective optimization for vibration reduction in composite plate structures using constrained layer damping. Computers & Structures, 2017: 105810

[2] Nie J, Wang J, Gou S, Zhu Y, Fan J. Technological development and engineering applications of novel steel-concrete composite structures. Frontiers of Structural and Civil Engineering, 2019, 13(1): 1–4

[3] Budarapu P R, Sudhir Sastry Y B, Natarajan R. Design concepts of an aircraft wing: Composite and morphing airfoil with auxetic structures. Frontiers of Structural and Civil Engineering, 2016, 10(4): 394–408

[4] Sakata S, Okuda K, Ikeda K. Stochastic analysis of laminated composite plate considering stochastic homogenization problem. Frontiers of Structural and Civil Engineering, 2015, 9(2): 141–153

[5] Asadi H, Kiani Y, Aghdam M M, Shakeri M. Enhanced thermal buckling of laminated composite cylindrical shells with shape memory alloy. Journal of Composite Materials, 2016, 50(2): 243–256

[6] Vo-Duy T, Ho-Huu V, Nguyen-Thoi T. Free vibration analysis of laminated FG-CNT reinforced composite beams using finite element method. Frontiers of Structural and Civil Engineering, 2018, 13(2): 324–336

[7] Fattahi I, Mirdamadi H R. A novel 3D skeletal frame topology for energy harvesting systems. Microelectronics Journal, 2019, 83: 6–17

[8] Rahmati M, Mirdamadi H R, Goli S. Divergence instability of pipes conveying fluid with uncertain flow velocity. Physica A, 2018, 491: 650–665

[9] Fattahi I, Mirdamadi H R. Novel composite finite element model for piezoelectric energy harvesters based on 3D beam kinematics. Composite Structures, 2017, 179: 161–171

[10] Fattahi I, Mirdamadi H R. A novel multimodal and multidirectional energy harvester by asymmetric 3D skeletal frame structures. Journal of the Brazilian Society of Mechanical Sciences and Engineering, 2020, 42(6): 274

[11] Naeeni IP, Keshavarzi A, Fattahi I. Parametric study on the geometric and kinetic aspects of the slider-crank mechanism. Iranian Journal of Science and Technology, Transactions of Mechanical Engineering, 2019, 43(3): 405–417

[12] Fattahi I, Mirdamadi H R. Electro-vibration modeling and response of 3D skeletal frame configuration for energy harvesters. Extreme Mechanics Letters, 2020: 100646

[13] Rabczuk T, Ren H, Zhuang X. A nonlocal operator method for partial differential equations with application to electromagnetic waveguide problem. Computers, Materials & Continua, 2019, 59(1): 31–55

[14] Anitescu C, Atroshchenko E, Alajlan N, Rabczuk T. Artificial neural network methods for the solution of second order boundary value problems. Computers, Materials & Continua, 2019, 59(1): 345–359

[15] Guo H, Zhuang X, Rabczuk T. A deep collocation method for the bending analysis of Kirchhoff plate. Computers Materials & Continua., 2019, 59(2): 433–456

[16] Hamdia K M, Ghasemi H, Zhuang X, Alajlan N, Rabczuk T. Sensitivity and uncertainty analysis for flexoelectric nanostructures. Computer Methods in Applied Mechanics and Engineering, 2018, 337: 95–109

[17] Ghasemi H, Park H S, Rabczuk T. A multi-material level set-based topology optimization of flexoelectric composites. Computer Methods in Applied Mechanics and Engineering, 2018, 332: 47–62

[18] Ghasemi H, Park H S, Rabczuk T. A level-set based IGA formulation for topology optimization of flexoelectric materials. Computer Methods in Applied Mechanics and Engineering, 2017, 313: 239–258

[19] Nanthakumar S S, Lahmer T, Zhuang X, Zi G, Rabczuk T. Detection of material interfaces using a regularized level set method in piezoelectric structures. Inverse Problems in Science and Engineering, 2016, 24(1): 153–176

[20] Grundberg S, Girhammar U A. A comparison of exact and approximate analyses of partially interacting composite beam-columns. In: The Fifteenth International Conference on Civil, Structural and Environmental Engineering Computing. Civil-Comp Press, 2015

[21] Wang Q, Shi D, Liang Q. Free vibration analysis of axially loaded laminated composite beams with general boundary conditions by using a modified Fourier-Ritz approach. Journal of Composite Materials, 2016, 50(15): 2111–2135

[22] Wenzel C, Vidal P, D’Ottavio M, Polit O. Coupling of heterogeneous kinematics and Finite Element approximations applied to composite beam structures. Composite Structures, 2014, 116: 177–192

[23] Polit O, Gallimard L, Vidal P, D’Ottavio M, Giunta G, Belouettar S. An analysis of composite beams by means of hierarchical finite elements and a variables separation method. Computers & Structures, 2015, 158: 15–29

[24] Karami M A, Inman D J. Analytical modeling and experimental verification of the vibrations of the zigzag microstructure for energy harvesting. Journal of Vibration and Acoustics, 2011, 133(1): 011002

[25] Ananda Babu A, Edwin Sudhagar P, Rajamohan V. Dynamic characterization of thickness tapered laminated composite plates. Journal of Vibration and Control, 2016, 22(16): 3555–3575

[26] Milazzo A. Refined equivalent single layer formulations and finite elements for smart laminates free vibrations. Composites. Part B, Engineering, 2014, 61: 238–253

[27] Panda S, Kumar A. A design of active constrained layer damping treatment for vibration control of circular cylindrical shell structure. Journal of Vibration and Control, 2016, 24: 1077546316670071

[28] Tornabene F. General higher-order layer-wise theory for free vibrations of doubly-curved laminated composite shells and panels. Mechanics of Advanced Materials and Structures, 2016, 23(9): 1046–1067

[29] Reddy J N. Mechanics of Laminated Composite Plates: Theory and Analysis. CRC press, 1997

[30] Weaver J R, Johnston P R. Structural Dynamics by Finite Elements. Englewood Cliffs, NJ: Prentice-Hall, 1987

[31] Petyt M. Introduction to Finite Element Vibration Analysis. Cambridge University Press, 2010

[32] Iannacci J, Sordo G. Up-scaled macro-device implementation of a MEMS wideband vibration piezoelectric energy harvester design concept. Microsystem Technologies, 2016, 22(7): 1639–1651

[33] Muthalif A G, Nordin N D. Optimal piezoelectric beam shape for single and broadband vibration energy harvesting: Modeling, simulation and experimental results. Mechanical Systems and Signal Processing, 2015, 54–55: 417–426

[34] Leo D J. Engineering Analysis of Smart Material Systems. Hoboken, NJ: John Wiley & Sons, 2007

[35] Erturk A, Inman D J. Piezoelectric Energy Harvesting. John Wiley & Sons, 2011

[36] Cook R D. Concepts and Applications of Finite Element Analysis. John Wiley & Sons, 2007

[37] Zhang Y, Zhuang X. Cracking elements method for dynamic brittle fracture. Theoretical and Applied Fracture Mechanics, 2019, 102: 1–9

[38] Zhang Y, Zhuang X. Cracking elements: A self-propagating strong discontinuity embedded approach for quasi-brittle fracture. Finite Elements in Analysis and Design, 2018, 144: 84–100

[39] Zhuang X, Wang Q, Zhu H. Multiscale modelling of hydro-mechanical couplings in quasi-brittle materials. International Journal of Fracture, 2017, 204(1): 1–27

[40] Ren H L, Zhuang X Y, Anitescu C, Rabczuk T. An explicit phase field method for brittle dynamic fracture. Computers & Structures, 2019, 217: 45–56

RIGHTS & PERMISSIONS

Higher Education Press