Nonlinear Finite Element Analysis of Free and Forced Vibrations of Cellular Plates Having Triangularly Prismatic Metamaterial Cores: A Strain Gradient Plate Model Approach


 The main objective of this paper is to develop a theoretically and numerically reliable and efficient methodology based on combining a finite element method and a strain gradient shear deformation plate model accounting for the nonlinear free and forced vibrations of cellular plates having equitriangularly prismatic metamaterial cores. The proposed model based on the nonlinear finite element strain gradient elasticity is developed for the first time to provide a computationally efficient framework for the simulation of the underlying nonlinear dynamics of cellular plates with advanced microarchitectures. The corresponding governing equations follow Mindlin’s SG elasticity theory including the micro-inertia effect applied to the first-order shear deformation plate theory along with the nonlinear von Kármán kinematics. Standard and higher-order computational homogenization methods determine the classical and strain gradient material constants, respectively. A higher-order \({C}^{1}\)-continuous 6-node finite element is adopted for the discretization of the governing variational formulation with respect to the spatial domain, and an arc-length continuation technique along with time periodic discretization is implemented to solve the resulting nonlinear time-dependent problem. Through a set of comparative studies with 3D full-field models as references, the accuracy and efficacy of the proposed dimension reduction methodology are demonstrated for a diverse range of problem parameters for analyzing the large-amplitude dynamic structural response.


Introduction
The rapid development of manufacturing technologies, in the field of additive manufacturing processes particularly, expands the spaces of production procedures and applications of advanced mechanical metamaterials and microarchitectural structures in various engineering fields [1]. According to these developments and the great fundamental characteristics of lattice architectures such as adjustable stiffness-to-weigh ratio and high surface-to-volume ratio, there are essential demands for developing reliable and efficient computational methodologies to predict the thermomechanical characteristics of these state-of-the-art artifacts.
The structural models based on the classical continuum mechanics have been effectively utilized to simulate the static and dynamic behavior of various engineering structures already for decades if not centuries. However, the complex internal microarchitectures of advanced cellular structures and the high computational costs of the common numerical approaches such as the conventional finite element methods, when adopted in the form of detailed modeling, considerably reduce the applicability of classical models regarding the structural simulation of mechanical metamaterials. In contrast, the generalized continuum theories like the microcontinuum theory [2], strain gradient theory (SGT) [3,4], and nonlocal theory [5] propose higher-order (or higher-grade) models by which the complexity of the microarchitecture-dependent behavior of cellular or lattice structures or the microstructure-dependent mechanics of micro/nanostructures can be efficiently handled. Indeed, the application of generalized continuum theories in the field of mechanics can be considered from two different viewpoints. First, by following experimental observations [6,7] on the size-dependent behavior of structures at micro/nano scales, these higher-order models have been extensively applied to capture this size-effect [8][9][10][11][12][13][14]. In the meantime, several studies on the linear and nonlinear structural models of beams [Error! Reference source not found. -17], plates [18-Error! Reference source not found.], and shells [23] within the SGT have been reported in the literature. Besides, different numerical solution methods such as higher-order finite element methods [24][25][26][27][28][29], and isogeometric analysis [30][31][32][33] have been developed. Second, the generalized continuum theories can be usefully employed to consider the microarchitecture-dependent mechanics of lattice or cellular structures [34][35][36][37][38][39][40]. Indeed, by using these higher-order theories, the internal microarchitecture can be homogenized via generalized constitutive relations, possibly together with dimension reduction models, which remarkably decreases both modeling and computation costs.
A literature survey discloses a rather limited number of studies on the mechanics of cellular materials under the higher-order continuum theories [34][35][36][37][38][39][40]. Since solids with different microarchitectures under various loading conditions and deformation states might not be appropriately modeled by using a unified generalized continuum theory, different models such as the Cosserat theory [41], micropolar theory [35,42,43] and higher gradient theories [34,[37][38][39][40]44,45,[46][47][48] have been developed. Closely connected to the present study, micropolar plate models were adopted in [42,43] to present a 2D equivalent singlelayer for lattice core sandwich plates and panels, and the resulting problem was solved by linear Lagrange interpolation finite elements along with the selective reduced integration approach. SGT-based theoretical studies on the structural analysis and the determination of higher-order constitutive parameters of cellular beams and plates with prismatic cores have been reported by Khakalo et al. in [46][47][48], with the corresponding kinematic formulations for both thin and thick structures and with isogeometric analysis for numerical solutions. Besides, the nonlinear static bending of cellular plates with an equitriangular microarchitecture was recently studied by Torabi and Niiranen in [49] highlighting the importance of geometric nonlinearity. Analytical studies on pantographic structures within a higher gradient theory with numerically and experimentally validated results can be found in [44,45]. Numerical homogenization within the SGT for a two-dimensional lattice structure was presented by Auffray et al. [38].
There exists only a few contributions on the dynamics of cellular structures within generalized continuum theories: the wave propagation analysis in [50,51]; the investigations of linear vibrations based on the micropolar theory [43,52,53], couple-stress theory [54], nonlocal theory [55] and SGT [46]. However, the effects of geometrical nonlinearities on the vibrations of cellular structures within the generalized continuum theories have not been studied. Accordingly, inspired by the work of Khakalo and Niiranen [48] in which the classical and SG constitutive parameters of SG plate models have been determined via numerical homogenization techniques, this study investigates by adopting a higher-order finite element (FE) method the geometrically nonlinear free/forced vibrations of cellular plates having triangularly prismatic cores formed by extruding 2D equitriangular lattices in the third coordinate direction.
This study formulates the problems of structural dynamics by following Mindlin's anisotropic SGT (Section 1) along with the first-order shear deformation plate theory (FSDPT) and von Kármán's hypothesis (Section 2). The bending-related classical and SG constitutive parameters from [48] are extended to the nonlinear formulation according to [49]. To facilitate the implementation of a higherorder FE method, a matrix form Hamiltonian of the problem set is first derived (Section 4), whereas the corresponding FE formulation (Section 5) is obtained by relying on a quasi-conforming 6-node triangular element (Section 6), to appropriately discretize the variational formulation with respect to the spatial domain. Time integration and arc-length continuation are finally employed for solving the nonlinear time-

Anisotropic strain gradient theory
On the basis of the fundamentals of the SG theory proposed by Mindlin [19], the strain energy of solids (with the occupied volume of ) is defined by using strain and first strain gradient tensors ( , ), and corresponding classical and higher-order stress tensors ( , Σ ) as where the kinematic and constitutive relations for the strain and classical Cauchy stress tensors are expressed by with as the components of the displacement and as the classical material stiffness tensor. In addition, the first strain gradient tensor and associated higher-order stress tensor are presented as with as the sixth-order stiffness tensor related to the SG theory. On the other hand, by considering as the inertial length scale parameters, the kinematic energy is expressed as
In accordance with the displacement field and by following von Karman's nonlinear kinematic relations, the strain vector is presented as We note that for instance 1,1 is the first-order derivative of 1 to coordinate 1 . With the aid of To make the governing equations more appropriate for a finite element study, the matrix forms of strain and strain gradient vectors are introduced by the following relations: . Besides, the nonlinear classical ( ) and strain gradient ( , = 1,2) matrix operators are expressed as ].
Furthermore, the matrices , and , (i=1,2) are given by the following relations: We note that 33 and 22 are the 3-by-3 and 2-by-2 identity matrices, respectively, 22 is the 2-by-2 zero matrix and 〈〉 signify the diag function.
On the other hand, by the employment of the Voigt notation for the constitutive laws defined in the second relations of Eqs. (2) and (3) where is the classical material stiffness matrix presented as with as the shear correction factor for the FSDPT. Besides, ( = 1,2,3,4) are the higher-order material stiffness matrices associated with the SG elasticity and can be given as follows: ].
The constants and in Eqs. (21) and (22) can be determined by using generalized numerical homogenization procedures for any specific microarchitecture lattices.

Hamilton's Principle
Hamilton's principle is applied to present the dynamic governing equations. By setting Π, , and , respectively, for strain energy, kinetic energy, and the work done by external loads, where stands for the first variation operator, it can be expressed in the form Under the strain energy of Eq. (1) and by following constitutive relations of Eq. (20), the total strain energy related to the SG theory are represented as It is worth noting that by considering the geometrically nonlinear kinematic relations of Eq. (1), the first variations of the strain and SG vectors are defined as follows: = ( + ) , Besides, by using Eq.
considering only the transverse loads = [0 0 0 0] indicating the vector of surface forces. The next section provides detailed formulations on the FE discretization procedure for the nonlinear vibration problem of the SG plate model.

Finite element discretization procedure
In accordance with the basic idea of the FE analysis, the displacement field can be approximated in each element by using appropriate approximate functions: where is the vector of the nodal unknowns related to the displacement components and where are the approximate functions and signifies the number of degrees of freedom (DOF) per node. By using the approximate expression presented for the displacement vector and by substituting Eq.
(30) into (12), the approximations of the strain and strain gradient vectors within each element are expressed by where the linear ( ) and nonlinear ( ) FE matrix operators of the classical theory are given as while for the SG counterparts, they can be written as Now, by substituting Eqs. (32) and (35) into the variation of the strain energy functional Eq. (24), FE approximations of the energy is presented as where the first and second integrals capture the classical and SG effects, respectively. The following material matrices, produced by integrations over the thickness direction, are written as where = ℕ.
Finally, substituting the fundamental equations given in Eqs. (36) and (38) into Hamilton's principle, results in the nonlinear FE governing equation: in which the stiffness matrix = + addresses both linear ( ) and nonlinear parts ( ) where the linear stiffness matrix consists of two terms, i.e., = + with where superscripts and , respectively, signify the classical and SG parts. In a similar manner, the nonlinear stiffness matrix can be divided into two parts, i.e., = + with In addition, the mass matrix and force vector are defined as By considering the harmonic excitation under a spatially uniform transverse load with an excitation frequency , and taking dissipation into account through the external viscous damping effect, the FE governing equation can be rewritten as follows: with as the damping matrix defined as where is the damping coefficient and is the linear fundamental frequency of the structure. The nonlinear FE governing equation provided in Eq. (44) can be solved by the numerical solution procedure described in [56,57] in which the time-periodic differential operators and an arc-length continuation technique are employed to present the nonlinear free/forced vibration responses. The stiffness and mass matrices as well as the force vector defined in Eqs. (40)- (43) are obtained by employing a quasiconforming six-node triangular element detailed next.

Quasi-conforming higher-order triangular element
In the majority of engineering applications, the standard 0 -continuous 2D and 3D elements can be effectively employed and, accordingly, different commercial FE software can be utilized for numerical simulations. However, SG elasticity requires, due to the existence of the higher-order derivations of the field variables in the energy functional, either 1 -continuous elements or 0 -continuous elements but a mixed formulation. The present study follows the former route by adopting a higher-order element: a quasi-conforming six-node triangular element.
To present the fundamentals of the element, the local area (LA) coordinate system of 1 , 2 , 3 is employed. Besides the value of the field variables, their first-order derivatives are considered as additional nodal values. To give more details about the element, let us consider as a scalar field. As explained above, the values of scalar field (Φ 0 ) and its first-order derivatives with respect to 1 and 2 (Φ 1 , Φ 2 ) are taken as the nodal values, as demonstrated in Figure 1. We note that stands for the number of nodes, and since a six-node triangle is considered, = 1,2, … ,6.
The associated mathematical formulations for the described nodal values are which can be represented in the vector form = ( ) .

Numerical results
The nonlinear shear deformation SG plate model, discretized by the quasi-conforming higher-order triangular finite element, will be next applied for studying the microarchitecture-dependent free and forced vibration of cellular plate structures. By employing the SG plate model, the mechanical behavior of three-dimensional cellular structures with prismatic lattice microarchitectures can be modeled by a two-dimensional model significantly reducing the computational cost of the numerical structural analysis.
An example of the cellular plates with an equitriangular lattice architecture is shown in Figure 2 and geometrical parameters are given in Table 1. The SG plate model with the constitutive parameters (as presented in the material matrices , 1 , 2 , 3 , 4 in Eqs. (20)- (22) and further in (24)) can be determined by using a homogenization technique in which the deformations of the 2D SG plate model are compared to the reference full-field model of the classical three-dimensional continuum mechanics. For the specific equitriangular lattice architecture considered in the present work, Kakalo and Niiranen [48] have determined the classical and SG material parameters for a reduced SG plate theory, augmented by Torabi and Niiranen [49]. The nonzero constitutive parameters are given in Tables 2 and 3. Since the cellular plate is assumed to be made of stainless steel, the mass density is = 7800 (kg/m 3 ), and by considering the geometrical parameters of the unit cell (see Figure 2b), the value of the factor used in Eq.      where the normal and tangential components of the displacements and rotation are given by the relations = 1 1 + 2 2 , = − 2 1 + 1 2 , Ψ = 1 1 + 2 2 and Ψ = − 2 1 + 1 2 with 1 and 2 as the unit normal vector components.
Before going through the details of the numerical results, the impacts of the micro inertia are studied to see how significantly this term affects the dynamic analysis of lattice structures. For this purpose, the variations of the linear frequency parameters (Ω = 2 /ℎ√ / 11 ) of the square lattice plate versus the mode numbers are presented in Figure 5 for different values of the intrinsic inertial length scale parameter γ (see Eqs. (4), (26) and (38)). The single-layer lattice plate with fully clamped boundary conditions is taken as an example. To give more quantitative results, the values of the first five frequency values are tabulated in Table 4. It is seen that the inertial length scale corresponding to the micro inertia term does not considerably affect the lower modes of the frequency parameter (at least for mode numbers lower than 10). Based on this finding, the micro inertia term is neglected in what follows, i.e., γ = 0.

Convergence analysis and model comparisons
This section provides some numerical results to express the accuracy and efficiency of the SG plate model

Microarchitecture effects on the nonlinear free/forced vibration responses
In this section, effects of the microarchitecture on the nonlinear free/forced vibration responses are studied. In this regard, variations of the non-dimensional deflection ( /ℎ) versus the frequency ratio are presented in Figures 18-20 as the nonlinear vibration responses for various shapes of clamped plates.
The non-dimensional force amplitude is defined as = 4 / 11 ℎ 4 . Different numbers of lattice layers are considered with a constant thickness ratio to highlight the importance of lattice architectures.
Comparing the free vibration responses for different numbers of lattice layers specifies that even with a similar thickness ratio, the one-layer plate experiences a stiffer dynamic behavior so that for the certain value of the non-dimensional deflection it has a smaller frequency ratio. Indeed, the plate with the onelayer lattice has a stiffer structure with a larger linear frequency so that it has a smaller frequency ratio.
Similar dynamic behavior can be observed for the forced vibration results where the one-layer plate has a lower peak in the frequency response following its stiffer behavior owing to the lattice microarchitecture effects. The results also reveal that by the increase of the lattice layers the importance of the microarchitecture decreases and for instance the nonlinear vibration responses of three-and four-layer plates are quite similar.

Conclusion
The microarchitecture-dependent large-amplitude free/forced vibrations of cellular plates with prismatic cores was studied within a strain gradient plate model discretized by finite element method. In the theoretical sections, the nonlinear equations of motion were derived by using the strain gradient elasticity theory and FSDPT along with the von Karman geometrically nonlinear kinematics. A detailed variational formulation was presented and a quasi-conforming higher-order triangular element was introduced to conduct the finite element discretization.
In the sections for numerical results, the efficiency of the finite element analysis was confirmed by convergence studies. Regarding the model, the effects of micro inertia were first then taken into account, but it was revealed that the inertial length scale parameter has no remarkable impact on the lowest natural frequencies investigated in the rest of the examples. For a wide range of examples of cellular plates with different shapes of mid-surface, the response of the plate model was compared to the response of the corresponding three-dimensional full-field model used as a reference, to address the accuracy of the proposed two-dimensional model. The results indicated that in the case of constant thickness ratio, the one-layer lattice plate has the stiffest mechanical behavior with higher natural frequencies due to the microarchitecture effects. It is also observed that for thin cellular plates the lattice microarchitecture has a more significant stiffening effect than for thick plates. Besides, it is found that due to the higher impact of the lattice microarchitecture, the one-layer plate has a lower peak of the forced vibration response.
Comparing the results for plates with different numbers of lattice layers revealed that by the increase of the number of layers the importance of the lattice microarchitecture diminishes and the dynamic behavior can be captured by the classical shear-deformable plate theory. All in all, the results provided on the linear and nonlinear vibration of cellular plates with equitriangular lattice core, along with our previously published results for plate bending statics [48,49] to be extended towards shell structuresprovide a comprehensive and reliable basis and efficient simulation tools for the structural design and engineering of microarchitectural plate structures.

Data Availability Statement:
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.