Two-Scale Asymptotic Homogenization Method for Composite Kirchhoff Plates with in-Plane Periodicity

Abstract

This paper develops a two-scale asymptotic homogenization method for periodic composite Kirchhoff plates. In this method, a three-dimensional (3D) periodic plate problem is simplified as a Kirchhoff plate problem, which is governed by a fourth-order uniformly elliptic partial differential equation (PDE) with periodically oscillating coefficients. Then, a two-scale solution in an asymptotic expansion form is presented for the PDE, and it is found that the first-order perturbed displacement in the asymptotic solution is zero. Additionally, periodic boundary and normalization constraint conditions are proposed to determine the unique solution to unit cell problems. Moreover, standard finite element formulations for calculating the perturbed displacements are derived from the principle of virtual work. Physical interpretations of the influence functions are presented by analyzing the properties of self-balanced quasi-load vectors used for solving the influence functions. Numerical comparisons show that the present method is physically acceptable and highly accurate.

1. Introduction

The development of modern technology is inseparable from composite materials. Composite material structures have been extensively used in various engineering fields in the past few decades. However, the heterogeneity of microstructures brings high computational costs to the modeling and simulation of composite structures. Since composite structures have various scale characteristics in nature, multiscale homogenization methods are an important tool for the analysis of composite structures from the point of view of accuracy and efficiency.

In recent decades, many multiscale homogenization methods have been proposed to deal with composite structures, such as the two-scale asymptotic homogenization method (AHM) [1,2,3], the multiscale eigenelement method (MEM) [4,5,6], the heterogeneous multiscale method (HMM) [7,8], the variational asymptotic method (VAM) [9,10], and for many other multiscale homogenization methods referred to [11,12] and the references cited therein. Among the above numerical homogenization methods [1,2,3,4,5,6,7,8,9,10], the AHM [1,2,3] is one of the most representative ones with a rigorous mathematical foundation and has been widely used in the homogenization analysis of periodic composite structures for statics [13,14,15,16,17,18] and dynamics [19,20,21,22]. However, the AHM [1,2,3,13,14,15,16,17,18,19,20,21,22] was mainly confined to coping with structural problems with omnidirectional periodicity, such as two-dimensional (2D) plane or three-dimensional (3D) solid problems, where the governing equations are second-order elliptic partial differential equations (PDEs) with periodically oscillating coefficients. Whereas, studies about the transverse homogenization of composite plate structures with only in-plane periodicity are still inadequate. A unique property of this kind of problem is that the plates deform in the thickness direction, which is perpendicular to the direction in which the periodicity exists. To realize the transverse homogenization of periodic plates, we have to choose one of the plate theories, which are briefly reviewed below, before reviewing the works related to the homogenization.

The Kirchhoff plate theory [23] is the simplest plate theory in which only deflection is the independent variable. Brunelle et al. [24] pointed out that the Kirchhoff plate theory is not suitable for moderately thick plates or those with a large ratio of elastic modulus to shear modulus, since the Kirchhoff plate theory ignores the transverse shear deformation. To improve the accuracy of the Kirchhoff plate theory, Reissner and Mindlin et al. [25,26] constructed the first-order shear plate theory by introducing two additional independent angles to describe the transverse shear deformations, but it cannot capture warping deformation, so the shear correction coefficient was introduced to update shear stiffness. Furthermore, several higher-order shear plate theories [27,28,29,30,31,32,33,34] were proposed by expanding in-plane displacements into higher-order power functions of thickness coordinates, one representative of which is the third-order shear deformation theory [32]. Additionally, layer-wise theories [35,36,37,38] provide a powerful tool to simulate the stress distribution of laminates, achieving higher accuracy compared with single-layer plate theories [39,40].

However, the existing plate theories [23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40] were developed for homogeneous plates or laminated composite plates. For composite plate structures with periodic microstructures, these plate theories cannot be used directly, and they should be used together with the solutions of the unit cell problem to effectively fulfill the static and dynamic analyses of the periodic plates.

To deal with periodic thin plate structures, some work focused on the simplification of 3D periodic thin plate problems to 2D homogenized Kirchhoff plate models in analytical manners. Libove et al. [41] were one of the earliest researchers to investigate the homogenized properties of sandwich plates and state that the sandwich panel can be transformed into an equivalent homogeneous panel with elastic constants. Briassoulis et al. [42] held that corrugated thin plates composed of anisotropic materials could be equivalent to homogenized orthotropic thin plates, and derived the analytical expressions of the effective bending stiffnesses in accordance with force-distortion relationships. Richard et al. [43], Samanta et al. [44], and Kress et al. [45] achieved the analytical results of effective stiffnesses of sinusoidal, trapezoidal, and circular corrugated plates. Based on the principle of strain energy equivalence, Xia et al. [46] established a homogenization-based analytical model to obtain the explicit expressions of effective in-plane tensile stiffness and out-of-plane bending stiffness for various sandwich plates. As for periodic moderately thick plate structures, transverse shear effects cannot be neglected, so the Reissner–Mindlin plate model should be considered in the process of stiffness predictions. According to the relation between force and deformation, Libove et al. [41] derived the effective bending stiffness and transverse shear stiffness of corrugated sandwich plates. Following the work of Libove et al. [41], Lok et al. [47,48], Fung et al. [49,50,51], and Leekitwattana et al. [52] derived the effective properties of some periodic, moderately thick plates in accordance with force-distortion relationships.

However, the abovementioned analytical approaches [41,42,43,44,45,46,47,48,49,50,51,52] to effective constant predictions may encounter difficulties when addressing periodic plates with complex microstructures. Hence, numerical homogenization methods may be a necessary choice in such situations. For composite plates with only in-plane periodic microstructures, Nasution et al. [53] used the AHM to achieve the in-plane effective properties by relieving the periodicity in the thickness direction. Nevertheless, the traditional AHM [1,2,3,13,14,15,16,17,18,19,20,21,22] cannot well capture out-of-plane deformations produced by transverse loads due to few or no repeated unit cells in the thickness direction. To achieve the out-of-plane effective stiffnesses of periodic thin plates, some work [54,55,56,57] simplified 3D plate problems to 2D homogenized Kirchhoff plate problems by assuming the behaviors in the thickness direction. However, these works [54,55,56,57] required the solution of complex 3D unit cell problems and were difficult to achieve the micro deformations under real loads.

Different from previous work, the main purpose of this work is to develop a homogenization method for composite thin plates with in-plane periodicity by combining the Kirchhoff plate theory with the two-scale AHM to predict the effective bending stiffness and achieve two-scale displacements. To the author’s knowledge, few studies have been conducted on the homogenization of the 2D periodic Kirchhoff plate. Kolpakov [58] applied the AHM to deal with the 2D inhomogeneous periodic plate with initial stresses. Recently, Faraci et al. [59] held that the two-scale AHM could be used to analyze sufficiently thin periodic plates based on the Kirchhoff plate theory, and the first two-order asymptotic solutions were presented. However, the above references [58,59] did not involve the solutions of higher perturbed terms or their physical interpretations.

In this work, we focus on the two-scale asymptotic solutions of 2D periodic Kirchhoff plates. The main novel contributions to this work are listed as follows:

(1)

The asymptotic expansion solutions of the fourth-order PDE with rapidly periodic oscillating coefficients were presented, and the accuracy of high-order perturbed terms was investigated quantitatively.

(2)

Constraint conditions for the unit cell problem were given and elaborated physically.

(3)

The influence functions were interpreted physically for a better understanding of the two-scale AHM.

(4)

The explicit analytically homogenized stiffness for periodic plates with layered structures was presented.

The rest of this paper is organized as follows: In Section 2, a 3D periodic plate problem is simplified as a 2D Kirchhoff plate model, and then the process of asymptotic homogenization for the 2D Kirchhoff plate model is presented; Section 3 provides the solutions for unit cell problems; in Section 4, the physical interpretations of influence functions are presented; in Section 5, numerical experiments are conducted to verify the proposed method; and conclusions are drawn in Section 6.

2. Asymptotic Homogenization of the Periodic Composite Plate

In this section, a 3D composite plate structure with in-plane periodicity is equivalent to a 2D composite plate according to the Kirchhoff plate theory, as shown in Figure 1. Then, the two-scale asymptotic homogenization method is employed to solve the 2D periodic Kirchhoff plate problem, i.e., a fourth-order elliptic PDE with periodically oscillating coefficients.

2.1. Dimensional Reduction to Composite Kirchhoff Plates with in-Plane Periodicity

According to the assumption of Kirchhoff plate theory, the displacement field of the plate for bending problems can be determined as:

$$\left\{\begin{array}{l}{u}_1^{\epsilon }=-x_3\frac{\partial w^{\epsilon }}{\partial x_1}\\ u_2^{\epsilon }=-x_3\frac{\partial w^{\epsilon }}{\partial x_2}\\ u_3^{\epsilon }=w^{\epsilon }\end{array}\right.$$

(1)

where $u_i^{\epsilon } \left(i=1, 2, 3\right)$ are the spatial displacements, and $w^{\epsilon }$ denotes the transverse displacement on the reference surface.

The elastic strain field can be evaluated from the displacement field as follows:

$${\mathit{\epsilon }}^{\epsilon }={\left[\begin{array}{ccc}{\epsilon }_1^{\epsilon }& {\epsilon }_2^{\epsilon }& {\gamma }_{12}^{\epsilon }\end{array}\right]}^{\mathrm{T}}={\left[\begin{array}{ccc}-x_3\frac{{\partial }^2{w}^{\epsilon }}{\partial x_1^2}& -x_3\frac{{\partial }^2{w}^{\epsilon }}{\partial x_2^2}& -2x_3\frac{{\partial }^2{w}^{\epsilon }}{\partial x_1\partial x_2}\end{array}\right]}^{\mathrm{T}}$$

(2)

where ${\epsilon }_1^{\epsilon }$ and ${\epsilon }_2^{\epsilon }$ denote the tensile strains along the ${\mathit{x}}_1$ and ${\mathit{x}}_2$ directions, respectively; ${\gamma }_{12}$ represents the shear strain of the plate.

The generalized strain vector ${\mathit{\kappa }}^{\epsilon }$ is defined as:

$${\mathbf{\kappa }}^{\epsilon }={\left[\begin{array}{ccc}{\kappa }_{11}^{\epsilon }& {\kappa }_{22}^{\epsilon }& 2{\kappa }_{12}^{\epsilon }\end{array}\right]}^{\mathrm{T}}={\left[\begin{array}{ccc}-\frac{{\partial }^2{w}^{\epsilon }}{\partial x_1^2}& -\frac{{\partial }^2{w}^{\epsilon }}{\partial x_2^2}& -2\frac{{\partial }^2{w}^{\epsilon }}{\partial x_1\partial x_2}\end{array}\right]}^{\mathrm{T}}$$

(3)

where ${\mathit{\kappa }}_{11}^{\epsilon }$ and ${\mathit{\kappa }}_{22}^{\epsilon }$ respectively represent the bending curvatures around the ${\mathit{x}}_2$ and ${\mathit{x}}_1$ axes, and ${\mathit{\kappa }}_{12}^{\epsilon }$ stands for the torsional curvature. From Equations (2) and (3), it shows that the elastic strain field depends on the generalized strains as follows:

$${\mathit{\epsilon }}^{\epsilon }=x_3{\mathbf{\kappa }}^{\epsilon }$$

(4)

The constitutive equation has the form as

$${\mathbf{\sigma }}^{\epsilon }={\mathit{E}}^{\epsilon }{\mathit{\epsilon }}^{\epsilon }$$

(5)

where

$${\mathbf{\sigma }}^{\epsilon }={\left[\begin{array}{ccc}{\sigma }_{11}^{\epsilon }& {\sigma }_{22}^{\epsilon }& {\sigma }_{12}^{\epsilon }\end{array}\right]}^{\mathrm{T}}$$

(6)

The definitions of internal moments are

$${\mathit{M}}^{\epsilon }={\displaystyle {\int }_{\partial h}{x}_3{\mathbf{\sigma }}^{\epsilon }\mathrm{d}{x}_3}={\mathit{D}}^{\epsilon }{\mathbf{\kappa }}^{\epsilon }$$

(7)

where

$${\mathit{M}}^{\epsilon }={\left[\begin{array}{ccc}{M}_{11}^{\epsilon }& M_{22}^{\epsilon }& M_{12}^{\epsilon }\end{array}\right]}^{\mathrm{T}}$$

(8)

$${\mathit{D}}^{\epsilon }={\displaystyle {\int }_{\partial h}{x}_3^2{\mathit{E}}^{\epsilon }\mathrm{d}{x}_3}$$

(9)

The equilibrium equations of an infinitesimal element of the plate with all the out-of-plane forces and externally applied loads, as shown in Figure 2, take the forms as:

$$\frac{\partial Q_{13}^{\epsilon }}{\partial x_1}+\frac{\partial Q_{23}^{\epsilon }}{\partial x_2}+q=0 \mathrm{in} \mathrm{\Omega }$$

(10)

$$\frac{\partial M_{11}^{\epsilon }}{\partial x_1}+\frac{\partial M_{21}^{\epsilon }}{\partial x_2}-Q_{13}^{\epsilon }=0 \mathrm{in} \mathrm{\Omega }$$

(11)

$$\frac{\partial M_{12}^{\epsilon }}{\partial x_1}+\frac{\partial M_{22}^{\epsilon }}{\partial x_2}-Q_{23}^{\epsilon }=0 \mathrm{in} \mathrm{\Omega }$$

(12)

where both $Q_{13}^{\epsilon }$ and $Q_{23}^{\epsilon }$ denote the transverse shear forces along the $x_3$, but they are defined in the cross sections perpendicular to the x₁ and x₂, respectively, as shown in Figure 2; $\mathrm{\Omega }$ denotes the domain of a 2D plate. Combined with Equations (10)–(12), the equilibrium equation of the Kirchhoff plate can be expressed as:

$$\frac{{\partial }^2{M}_{ij}^{\epsilon }}{\partial x_i\partial x_j}=-q \mathrm{in} \mathrm{\Omega }$$

(13)

where i, j as well as the following k and l assume values 1 and 2. Note that the summation convention is assumed for repeated indices throughout this paper. With Equation (7), the equilibrium Equation (13) can be rewritten as:

$$\frac{{\partial }^2}{\partial x_i\partial x_j}\left(D_{ijkl}^{\epsilon }\frac{{\partial }^2{w}^{\epsilon }}{\partial x_k\partial x_l}\right)=q \mathrm{in} \mathrm{\Omega }$$

(14)

Here, $D_{ijkl}^{\epsilon }$ denotes the bending stiffness tensor. For periodic Kirchhoff plates, $D_{ijkl}^{\epsilon }$ oscillates periodically in the plane, implying that the solution of Equation (14) is also oscillating. Hence, the following is to employ the two-scale asymptotic homogenization method to solve the fourth-order elliptic equation for periodic Kirchhoff plates.

2.2. Two-Scale Asymptotic Homogenization Method for 2D Periodic Kirchhoff Plates

Taking into account the small periodic configurations of Kirchhoff plates, this subsection first establishes a two-scale model with macroscopic variable $\mathit{x}$ and microscopic variable $\mathit{y}=\mathit{x}/\epsilon$. Here, the small parameter $\epsilon$ is the separation scale, and $\mathit{x}$ gives the position of a point in $\mathsf{\Omega }$, while $\mathit{y}$ denotes the position in a unit cell domain Y.

In the two-scale asymptotic homogenization method, the actual displacement ${\mathit{w}}^{\epsilon }$, or the solution of Equation (14), has an asymptotic expansion form with two scales as:

$$w^{\epsilon }(\mathit{x})=w_0(\mathit{x})+\epsilon w_1(\mathit{x},\mathit{y})+{\epsilon }^2{w}_2(\mathit{x},\mathit{y})+{\epsilon }^3{w}_3(\mathit{x},\mathit{y})+\cdots$$

(15)

where ${\mathit{w}}_0\left(\mathit{x}\right)$ denotes the homogenized displacement; ${\mathit{w}}_n(\mathit{x}, \mathit{y})$ (n = 1, 2, ⋯) represents the perturbed displacements caused by the heterogeneity of microstructures.

Let Ψ = Ψ($\mathit{x}, \mathit{y}$) be a function depending on two scales and then we have:

$${\Psi }^{\epsilon }\left(\mathit{x}\right)=\Psi \left(\mathit{x},\frac{\mathit{x}}{\epsilon }\right)$$

(16)

Note that

$$\frac{\partial {\Psi }^{\epsilon }}{\partial \mathit{x}}\left(\mathit{x}\right)=\frac{\partial \Psi }{\partial \mathit{x}}\left(\mathit{x},\frac{\mathit{x}}{\epsilon }\right)+\frac{1}{\epsilon }\frac{\partial \Psi }{\partial \mathit{y}}\left(\mathit{x},\frac{\mathit{x}}{\epsilon }\right)$$

(17)

$$\frac{{\partial }^2{\Psi }^{\epsilon }}{\partial {\mathit{x}}^2}\left(\mathit{x}\right)=\frac{{\partial }^2\Psi }{\partial {\mathit{x}}^2}\left(\mathit{x},\frac{\mathit{x}}{\epsilon }\right)+\frac{2}{\epsilon }\frac{{\partial }^2\Psi }{\partial \mathit{x}\partial \mathit{y}}\left(\mathit{x},\frac{\mathit{x}}{\epsilon }\right)+\frac{1}{{\epsilon }^2}\frac{{\partial }^2\Psi }{\partial {\mathit{y}}^2}\left(\mathit{x},\frac{\mathit{x}}{\epsilon }\right)$$

(18)

Substituting Equation (18) into Equation (14) results in the following:

$$\left(\frac{1}{{\epsilon }^4}{A}_0+\frac{2}{{\epsilon }^3}{A}_1+\frac{1}{{\epsilon }^2}{A}_2+\frac{2}{\epsilon }{A}_3+A_4\right)w^{\epsilon }\left(\mathit{x}\right)=q \mathrm{in} \mathrm{\Omega }$$

(19)

where

$$\left\{\begin{array}{l}{A}_0=\frac{{\partial }^2}{\partial y_i\partial y_j}\left(D_{ijkl}^{\epsilon }\frac{{\partial }^2}{\partial y_k\partial y_l}\right)\\ A_1=\frac{{\partial }^2}{\partial x_i\partial y_j}\left(D_{ijkl}^{\epsilon }\frac{{\partial }^2}{\partial y_k\partial y_l}\right)+\frac{{\partial }^2}{\partial y_i\partial y_j}\left(D_{ijkl}^{\epsilon }\frac{{\partial }^2}{\partial x_k\partial y_l}\right)\\ A_2=\frac{{\partial }^2}{\partial x_i\partial x_j}\left(D_{ijkl}^{\epsilon }\frac{{\partial }^2}{\partial y_k\partial y_l}\right)+\frac{{\partial }^2}{\partial y_i\partial y_j}\left(D_{ijkl}^{\epsilon }\frac{{\partial }^2}{\partial x_k\partial x_l}\right)+\frac{4{\partial }^2}{\partial x_i\partial y_j}\left(D_{ijkl}^{\epsilon }\frac{{\partial }^2}{\partial x_k\partial y_l}\right)\\ A_3=\frac{{\partial }^2}{\partial x_i\partial x_j}\left(D_{ijkl}^{\epsilon }\frac{{\partial }^2}{\partial x_k\partial y_l}\right)+\frac{{\partial }^2}{\partial x_i\partial y_j}\left(D_{ijkl}^{\epsilon }\frac{{\partial }^2}{\partial x_k\partial x_l}\right)\\ A_4=\frac{{\partial }^2}{\partial x_i\partial x_j}\left(D_{ijkl}^{\epsilon }\frac{{\partial }^2}{\partial x_k\partial x_l}\right)\end{array}\right.$$

(20)

By substituting Equation (15) into Equation (19) and equating the power-like terms of ε, a series of governing equations for unit cells are obtained as follows:

$$O\left({\epsilon }^{-4}\right): A_0{w}_0=0 \mathrm{in} \mathrm{\Omega }$$

(21)

$$O\left({\epsilon }^{-3}\right): A_0{w}_1=-2A_1{w}_0 \mathrm{in} \mathrm{\Omega }$$

(22)

$$O\left({\epsilon }^{-2}\right): A_0{w}_2=-2A_1{w}_1-A_2{w}_0 \mathrm{in} \mathrm{\Omega }$$

(23)

$$O\left({\epsilon }^{-1}\right): A_0{w}_3=-2A_1{w}_2-A_2{w}_1-2A_3{w}_0 \mathrm{in} \mathrm{\Omega }$$

(24)

$$O\left(1\right): A_0{w}_4=q-2A_1{w}_3-A_2{w}_2-2A_3{w}_1-A_4{w}_0 \mathrm{in} \mathrm{\Omega }$$

(25)

$$O\left({\epsilon }^n\right): A_0{w}_{n+4}=-2A_1{w}_{n+3}-A_2{w}_{n+2}-2A_3{w}_{n+1}-A_4{w}_n, (n\ge 1) \mathrm{in} \mathrm{\Omega }$$

(26)

from which the unknowns ${\mathit{w}}_n$ can be determined successively. The spatial periodicity of material configuration indicates the same periodicity for ${\mathit{w}}_n(\mathit{x}, \mathit{y})$ (n = 0, 1, 2, ⋯) in Y as follows:

$$w_n(\mathit{x}, \mathit{y})=w_n(\mathit{x}, \mathit{y}+\mathit{Y})$$

(27)

With the in-plane periodicity of $D_{ijkl}^{\epsilon }$ in $\mathit{y}$ and Equation (27), one can find from Equation (21) that ${\mathit{w}}_0$ is independent of $\mathit{y}$, showing the rationality of taking ${\mathit{w}}_0$ as the homogenized term in Equation (15), that is ${\mathit{w}}_0={\mathit{w}}_0\left(\mathit{x}\right)$. Then Equation (22) can be simplified into the same form as Equation (21), as

$$A_0{w}_1=0 \mathrm{in} \mathrm{\Omega }$$

(28)

It follows from Equation (28) that ${\mathit{w}}_1$ can only be a linear function of $\mathit{y}$. However, $D_{ijkl}^{\epsilon }$ and ${\mathit{w}}_1$ are periodic in $\mathit{y}$, so one can infer that ${\mathit{w}}_1$ must be independent of $\mathit{y}$. Note that in the asymptotic expansion expression (15), ${\mathit{w}}_0\left(\mathit{x}\right)$ represents the homogenized deformation, while ${\mathit{w}}_n(\mathit{x}, \mathit{y})$ (n = 1, 2, ⋯) denotes the perturbation displacements. Thus, we have the following:

$$w_1\left(\mathit{x}, \mathit{y}\right)=0$$

(29)

By substituting Equation (29) and ${\mathit{w}}_0$ into Equation (23), we have the governing equation concerning ${\mathit{w}}_2$ as:

$$\frac{{\partial }^2}{\partial y_i\partial y_j}\left(D_{ijkl}^{\epsilon }\frac{{\partial }^2{w}_2}{\partial y_k\partial y_l}\right)=-\frac{{\partial }^2}{\partial y_i\partial y_j}\left(D_{ijkl}^{\epsilon }\frac{{\partial }^2{w}_0}{\partial x_k\partial x_l}\right) \mathrm{in} \mathrm{\Omega }$$

(30)

The separation-of-variable method is employed to ${\mathit{w}}_2$ as follows:

$$w_2(\mathit{x}, \mathit{y})={\chi }_2^{mn}(\mathit{y})\frac{{\partial }^2{w}_0}{\partial x_m\partial x_n}$$

(31)

where ${\mathit{\chi }}_2^{mn}$ denotes the second-order influence function with symmetry on indices m and n, and we denote ${\mathit{\chi }}_2=\left[{\mathit{\chi }}_2^{11} {\mathit{\chi }}_2^{22} {\mathit{\chi }}_2^{12}\right]$ following the symmetry.

By substituting Equation (31) into Equation (30), the governing equation for ${\mathit{\chi }}_2^{mn}$ is obtained as follows:

$$\frac{{\partial }^2}{\partial y_i\partial y_j}\left(D_{ijkl}^{\epsilon }\frac{{\partial }^2{\chi }_2^{mn}}{\partial y_k\partial y_l}\right)=-\frac{{\partial }^2{D}_{ijmn}^{\epsilon }}{\partial y_i\partial y_j} \mathrm{in} Y$$

(32)

Similarly, ${\mathit{w}}_3$ in Equation (24) can also be expressed in a separation-of-variable form as follows:

$$w_3(\mathit{x}, \mathit{y})={\chi }_3^{mnp}(\mathit{y})\frac{{\partial }^3{w}_0}{\partial x_m\partial x_n\partial x_p}$$

(33)

where the third-order influence function ${\mathit{\chi }}_3^{mnp}$ is governed by:

$$\frac{{\partial }^2}{\partial y_i\partial y_j}\left(D_{ijkl}^{\epsilon }\frac{{\partial }^2{\chi }_3^{mnp}}{\partial y_k\partial y_l}\right)=-\frac{2\partial D_{ipmn}^{\epsilon }}{\partial y_i}-\frac{2\partial }{\partial y_i}\left(D_{ipkl}^{\epsilon }\frac{{\partial }^2{\chi }_2^{mn}}{\partial y_k\partial y_l}\right) -\frac{2{\partial }^2}{\partial y_i\partial y_j}\left(D_{ijkp}^{\epsilon }\frac{\partial {\chi }_2^{mn}}{\partial y_k}\right) \mathrm{in} Y$$

(34)

which is obtained by combining Equations (24), (29), (31), and (33) and ${\mathit{w}}_0.$ Note that ${\mathit{\chi }}_3^{mnp}$ also has a symmetry on indices m and n, i.e., ${\mathit{\chi }}_3=\left[{\mathit{\chi }}_3^{111} {\mathit{\chi }}_3^{221} {\mathit{\chi }}_3^{121} {\mathit{\chi }}_3^{112} {\mathit{\chi }}_3^{222} {\mathit{\chi }}_3^{122}\right]$ following the symmetry.

By integrating Equation (25) over the domain Y, one has

$$\frac{1}{\left|Y\right|}{\displaystyle {\int }_Y{A}_0{w}_4\mathrm{d}Y}=\frac{1}{\left|Y\right|}{\displaystyle {\int }_Y\left(q-2A_1{w}_3-A_2{w}_2-2A_3{w}_1-A_4{w}_0\right)\mathrm{d}Y}$$

(35)

where $\left|Y\right|$ is the area of a unit cell of a 2D plate. Following the periodicity of ${\mathit{w}}_n$ in Equation (27), and with Equations (29) and (31) together, Equation (35) can be simplified as follows:

$$\frac{{\partial }^2}{\partial x_i\partial x_j}\left(D_{ijkl}^{\mathrm{H}}\frac{{\partial }^2{w}_0}{\partial x_k\partial x_l}\right)=q \mathrm{in} \mathrm{\Omega }$$

(36)

where the homogenized stiffness $D_{ijkl}^{\mathrm{H}}$ is as shown:

$$D_{ijkl}^{\mathrm{H}}=\frac{1}{\left|Y\right|}{\displaystyle {\int }_Y\left(D_{ijkl}^{\epsilon }+D_{ijmn}^{\epsilon }\frac{{\partial }^2{\chi }_2^{kl}}{\partial y_m\partial y_n}\right)dY}$$

(37)

It follows from Equation (37) that $D_{ijkl}^{\mathrm{H}}$ depends on the second-order derivative of ${\mathit{\chi }}_2^{kl}$ and also the material distributions within the unit cell. The ${\mathit{\chi }}_2^{kl}$ can be obtained by solving the unit cell problems with the proposed constraint conditions in Section 3. After that, the homogenized solutions can be achieved with the obtained homogenized stiffness. Note that the well-posedness and convergence of the homogenized solution to the fourth-order PDE can be proved by the Lax–Milgram theorem and the two-scale convergence method, refer to [60,61,62,63].

3. Solutions for Unit Cell Problems

In this section, constraint conditions are first proposed to determine influence functions which are used to calculate the homogenized stiffness in Equation (37) and the two-scale solution in Equation (15); then, the analytical homogenized stiffnesses are presented for periodic layered structures; finally, finite element formulations are provided for unit cell problems.

3.1. Constraint Conditions for Unit Cell Problems

A key factor for obtaining the influence functions ${\mathit{\chi }}_n$ (n = 2, 3) is to determine the constraint conditions for unit cell problems (32) and (34). Since the perturbed term ${\mathit{w}}_n(\mathit{x}, \mathit{y})$ (n = 1, 2, ⋯) is periodic in Y, implying that ${\mathit{\chi }}_n$ (n = 2, 3) has the same periodicity in y as follows:

$$\left\{\begin{array}{l}{\left.{\mathit{\chi }}_n\right|}_{\partial Y_{1-}}={\left.{\mathit{\chi }}_n\right|}_{\partial Y_{1+}}\\ {\left.{\mathit{\chi }}_n\right|}_{\partial Y_{2-}}={\left.{\mathit{\chi }}_n\right|}_{\partial Y_{2+}}\end{array}\right.$$

(38)

where $\partial Y_{1\pm }$ and $\partial Y_{2\pm }$ respectively represent two pairs of periodic boundaries, as shown in Figure 3:

Equation (38) implies the continuity conditions of the interface displacements. Likewise, to satisfy the continuity conditions of normal rotation angles, normal bending moments, as well as equivalent shear forces, another three periodic conditions concerning ${\mathit{\chi }}_n$ (n = 2, 3) is required as follows:

$$\left\{\begin{array}{l}{\left.\frac{\partial {\mathit{\chi }}_n}{\partial y_1}\right|}_{\partial Y_{1-}}={\left.\frac{\partial {\mathit{\chi }}_n}{\partial y_1}\right|}_{\partial Y_{1+}}\\ {\left.\frac{\partial {\mathit{\chi }}_n}{\partial y_2}\right|}_{\partial Y_{2-}}={\left.\frac{\partial {\mathit{\chi }}_n}{\partial y_2}\right|}_{\partial Y_{2+}}\end{array}\right.$$

(39)

$$\left\{\begin{array}{l}{\left.{\mathit{M}}_{11}^n\right|}_{\partial Y_{1-}}={\left.{\mathit{M}}_{11}^n\right|}_{\partial Y_{1+}}, {\mathit{M}}_{11}^n=-D_{11kl}^{\epsilon }\frac{{\partial }^2{\mathit{\chi }}_n}{\partial y_k\partial y_l}\\ {\left.{\mathit{M}}_{22}^n\right|}_{\partial Y_{2-}}={\left.{\mathit{M}}_{22}^n\right|}_{\partial Y_{2+}}, {\mathit{M}}_{22}^n=-D_{22kl}^{\epsilon }\frac{{\partial }^2{\mathit{\chi }}_n}{\partial y_k\partial y_l}\end{array}\right.$$

(40)

$$\left\{\begin{array}{l}{\left.{\mathit{V}}_{13}^n\right|}_{\partial Y_{1-}}={\left.{\mathit{V}}_{13}^n\right|}_{\partial Y_{1+}}, {\mathit{V}}_{13}^n={\mathit{Q}}_{13}^n+\frac{\partial {\mathit{M}}_{12}^n}{\partial y_2}=-\frac{\partial }{\partial y_1}\left(D_{11kl}^{\epsilon }\frac{{\partial }^2{\mathit{\chi }}_n}{\partial y_k\partial y_l}\right)-\frac{2\partial }{\partial y_2}\left(D_{12kl}^{\epsilon }\frac{{\partial }^2{\mathit{\chi }}_n}{\partial y_k\partial y_l}\right)\\ {\left.{\mathit{V}}_{23}^n\right|}_{\partial Y_{2-}}={\left.{\mathit{V}}_{23}^n\right|}_{\partial Y_{2+}}, V_{23}^n={\mathit{Q}}_{23}^n+\frac{\partial {\mathit{M}}_{12}^n}{\partial y_1}=-\frac{\partial }{\partial y_2}\left(D_{22kl}^{\epsilon }\frac{{\partial }^2{\mathit{\chi }}_n}{\partial y_k\partial y_l}\right)-\frac{2\partial }{\partial y_1}\left(D_{12kl}^{\epsilon }\frac{{\partial }^2{\mathit{\chi }}_n}{\partial y_k\partial y_l}\right)\end{array}\right.$$

(41)

It is worth noting that periodic boundary conditions (38) eliminate the rigid rotation modes with respect to the $y_2$ and $y_1$ axes. In order to remove the remaining rigid translation mode along the transverse direction, ${\mathit{\chi }}_n$ (n = 2, 3) is required to satisfy the following normalization constraint condition:

$$\frac{1}{\left|Y\right|}{\displaystyle {\int }_Y{\mathbf{\chi }}_{n}dY}=\mathbf{0}$$

(42)

According to the Kirchhoff plate theory, C¹ continuity is supposed to be satisfied at the interface of unit cells for the periodic plate, and thus, another two periodic conditions are used to ensure the continuity of tangential rotation angles as follows:

$$\left\{\begin{array}{l}{\left.\frac{\partial {\mathit{\chi }}_n}{\partial y_2}\right|}_{\partial Y_{1-}}={\left.\frac{\partial {\mathit{\chi }}_n}{\partial y_2}\right|}_{\partial Y_{1+}}\\ {\left.\frac{\partial {\mathit{\chi }}_n}{\partial y_1}\right|}_{\partial Y_{2-}}={\left.\frac{\partial {\mathit{\chi }}_n}{\partial y_1}\right|}_{\partial Y_{2+}}\end{array}\right.$$

(43)

With the proposed constraint conditions (38)–(43), one can solve the unit cell problems (32) and (34) by analytical methods for simple structures or by the finite element method (FEM) for general structures. It is worth noting that the perturbed displacements ${\mathit{w}}_n$ (n = 2, 3) are defined as the products of the influence function ${\mathit{\chi }}_n$ and the corresponding n-th derivative of ${\mathit{w}}_0$. Here, ${\mathit{\chi }}_n$ reflecting the inhomogeneity of unit cells, is the solution of the unit cell problems (32) and (34), while ${\mathit{w}}_0$ is the homogenized (macroscopic) solution to the global problem (36) depending on the obtained homogenized stiffness $D_{ijkl}^{\mathrm{H}}$ and the external load q. As a result, only ${\mathit{w}}_0$ meets the boundary conditions of the plate, whereas ${\mathit{w}}_n$ cannot in general due to the use of periodic boundary conditions.

3.2. Analytical Homogenized Stiffnesses for Periodic Layered Plates

This section presents the analytical solution of homogenized stiffness for periodic layered plates, where the material properties change only in the $y_1$ direction, refer to Figure 4b as an example. Hence, $D_{ijkl}^{\epsilon }$ is the function of $y_1$, that is

$$D_{ijkl}^{\epsilon }\left(y_1,y_2\right)=D_{ijkl}^{\epsilon }\left(y_1\right) \mathrm{in} Y$$

(44)

Taking into account the constraint conditions (38)–(43) and the distribution of $D_{ijkl}^{\epsilon }$ as Equation (44), the solution of the cell problem (32) is assumed as a function of $y_1$, and has the following form:

$${\mathit{\chi }}_2={\mathit{\chi }}_2\left(y_1\right) \mathrm{in} Y$$

(45)

By substituting Equation (45) into (32), the unit cell problem (32) becomes a fourth-order ordinary differential equation as follows:

$$\frac{{\mathrm{d}}^2}{\mathrm{d}{y}_1^2}\left(D_{1111}^{\epsilon }\frac{{\mathrm{d}}^2{\chi }_2^{mn}}{\mathrm{d}{y}_1^2}\right)=-\frac{{\mathrm{d}}^2{D}_{11mn}^{\epsilon }}{\mathrm{d}{y}_1^2} \mathrm{in} Y$$

(46)

Inserting Equation (45) into Equation (37) results in:

$$D_{ijkl}^{\mathrm{H}}=\frac{1}{\left|Y\right|}{\displaystyle {\int }_Y\left(D_{ijkl}^{\epsilon }+D_{ij11}^{\epsilon }\frac{{\mathrm{d}}^2{\chi }_2^{kl}}{\mathrm{d}{y}_1^2}\right)dY}$$

(47)

Thus, the homogenized stiffness can be achieved if ${\mathrm{d}}^2{\chi }_2^{kl}/\mathrm{d}{y}_1^2$ is determined. With the periodic conditions (39) and (40), ${\mathrm{d}}^2{\chi }_2^{kl}/\mathrm{d}{y}_1^2$ can be obtained from Equation (46), as

$$\frac{{\mathrm{d}}^2{\chi }_2^{mn}}{\mathrm{d}{y}_1^2}=\frac{1}{D_{1111}^{\epsilon }}\left(\frac{〈D_{11mn}^{\epsilon }/D_{1111}^{\epsilon }〉}{〈D_{1111}^{\epsilon }{}^{-1}〉}-D_{11mn}^{\epsilon }\right) \mathrm{in} Y$$

(48)

where $\langle ·\rangle =\frac{1}{\left|Y\right|}{{\displaystyle \int }}_Y^{ }·\mathrm{d}Y$.

By substituting Equation (48) into Equation (47), one can obtain the explicit expressions of $D_{ijkl}^{\mathrm{H}}$ as:

$$\left\{\begin{array}{l}{D}_{1111}^H=〈D_{1111}^{\epsilon }{}^{-1}〉{}^{-1}\\ D_{1122}^H=〈\frac{D_{1122}^{\epsilon }}{D_{1111}^{\epsilon }}〉{〈D_{1111}^{\epsilon }{}^{-1}〉}^{-1}\\ D_{1112}^H=〈\frac{D_{1112}^{\epsilon }}{D_{1111}^{\epsilon }}〉{〈D_{1111}^{\epsilon }{}^{-1}〉}^{-1}\\ D_{2222}^H=〈\frac{D_{2211}^{\epsilon }}{D_{1111}^{\epsilon }}〉〈\frac{D_{1122}^{\epsilon }}{D_{1111}^{\epsilon }}〉{〈D_{1111}^{\epsilon }{}^{-1}〉}^{-1}+〈D_{2222}^{\epsilon }-\frac{D_{2211}^{\epsilon }{D}_{1122}^{\epsilon }}{D_{1111}^{\epsilon }}〉\\ D_{2212}^H=〈\frac{D_{2211}^{\epsilon }}{D_{1111}^{\epsilon }}〉〈\frac{D_{1112}^{\epsilon }}{D_{1111}^{\epsilon }}〉{〈D_{1111}^{\epsilon }{}^{-1}〉}^{-1}+〈D_{2212}^{\epsilon }-\frac{D_{2211}^{\epsilon }{D}_{1112}^{\epsilon }}{D_{1111}^{\epsilon }}〉\\ D_{1212}^H=〈\frac{D_{1211}^{\epsilon }}{D_{1111}^{\epsilon }}〉〈\frac{D_{1112}^{\epsilon }}{D_{1111}^{\epsilon }}〉{〈D_{1111}^{\epsilon }{}^{-1}〉}^{-1}+〈D_{1212}^{\epsilon }-\frac{D_{1211}^{\epsilon }{D}_{1112}^{\epsilon }}{D_{1111}^{\epsilon }}〉\end{array}\right.$$

(49)

It can be seen from Equation (49) that $D_{ijkl}^{\mathrm{H}}$ depends on the proportions of materials of layered structures in a highly nonlinear way.

3.3. Finite Element Formulations for Solving Unit Cell Problems

This section aims to provide the finite element formulations of the boundary value problems for ${\mathit{\chi }}_2$ and ${\mathit{\chi }}_3$. The Lagrange multiplier method [64] was used to consider the normalization constraint condition (42), and the master-slave method [65] was employed to deal with the periodic boundary conditions (38), (39) and (43).

It follows from the governing equations of ${\mathit{\chi }}_2$ and ${\mathit{\chi }}_3$ in Equations (32) and (34) that the right-side terms are actually self-balanced quasi-loads depending on the material constants. According to the principle of virtual work and the Lagrange multiplier method, the unit cell problems (32) and (34) with the normalization constraint condition (42) can be respectively written as:

$${\displaystyle {\int }_Y{D}_{ijkl}^{\epsilon }\frac{{\partial }^2{\chi }_2^{mn}}{\partial y_i\partial y_j}\frac{{\partial }^2\mathsf{\delta }{\chi }_2^{}}{\partial y_k\partial y_l}\mathrm{d}Y+}{\displaystyle {\int }_Y{D}_{ijmn}^{\epsilon }\frac{{\partial }^2\mathsf{\delta }{\chi }_2^{}}{\partial y_i\partial y_j}\mathrm{d}Y}+{\lambda }_2^{mn}{\displaystyle {\int }_Y\mathsf{\delta }{\chi }_2^{}\mathrm{d}Y}+\mathsf{\delta }{\lambda }_2{\displaystyle {\int }_Y{\chi }_2^{mn}\mathrm{d}Y}=0$$

(50)

$$\begin{array}{l}{\displaystyle {\int }_Y{D}_{ijkl}^{\epsilon }\frac{{\partial }^2{\chi }_3^{mnp}}{\partial y_i\partial y_j}\frac{{\partial }^2\mathsf{\delta }{\chi }_3^{}}{\partial y_k\partial y_l}\mathrm{d}Y}+{\displaystyle {\int }_Y\left(2D_{ipmn}^{\epsilon }\frac{\partial \mathsf{\delta }{\chi }_3^{}}{\partial y_i}-2D_{ipkl}^{\epsilon }\frac{\partial \mathsf{\delta }{\chi }_3^{}}{\partial y_i}\frac{{\partial }^2{\chi }_2^{mn}}{\partial y_k\partial y_l}+2D_{ijkp}^{\epsilon }\frac{{\partial }^2\mathsf{\delta }{\chi }_3^{}}{\partial y_i\partial y_j}\frac{\partial {\chi }_2^{mn}}{\partial y_k}\right)\mathrm{d}Y}\\ +{\lambda }_3^{mnp}{\displaystyle {\int }_Y\mathsf{\delta }{\chi }_3^{}\mathrm{d}Y}+\mathsf{\delta }{\lambda }_3{\displaystyle {\int }_Y{\chi }_3^{mnp}\mathrm{d}Y}=0\end{array}$$

(51)

where ${\mathit{\lambda }}_2^{mn}$ and ${\mathit{\lambda }}_3^{mnp}$ denote the Lagrange multipliers.

Considering the symmetric property of $D_{ijkl}^{\epsilon }$ and the arbitrariness of the virtual displacements $\mathsf{\delta }{\chi }_2$ and $\mathsf{\delta }{\chi }_3$, Equations (50) and (51) can be respectively discretized in matrix form as:

$$\left[\begin{array}{cc}\mathit{K}& {\mathit{C}}^{\mathrm{T}}\\ \mathit{C}& 0\end{array}\right]\left[\begin{array}{c}{\mathit{{\rm X}}}_2\\ {\mathit{\lambda }}_2\end{array}\right]=\left[\begin{array}{c}{\mathit{F}}_2\\ \mathbf{0}\end{array}\right]$$

(52)

$$\left[\begin{array}{cc}\mathit{K}& {\mathit{C}}^{\mathrm{T}}\\ \mathit{C}& 0\end{array}\right]\left[\begin{array}{c}{\mathit{{\rm X}}}_3\\ {\mathit{\lambda }}_3\end{array}\right]=\left[\begin{array}{c}{\mathit{F}}_3\\ \mathbf{0}\end{array}\right]$$

(53)

where

$$\mathit{K}={\displaystyle \sum _{e=1}{\displaystyle {\int }_{Y^e}{\mathit{B}}_2^T{\mathit{D}}^{\epsilon }{\mathit{B}}_2\mathrm{d}Y}}$$

(54)

$$\mathit{C}={\displaystyle \sum _{e=1}{\displaystyle {\int }_{Y^e}\mathit{N}\mathrm{d}Y}}$$

(55)

$${\mathit{{\rm X}}}_2={\displaystyle \sum _{e=1}{\mathit{{\rm X}}}_2^e}, {\mathit{{\rm X}}}_3={\displaystyle \sum _{e=1}{\mathit{{\rm X}}}_3^e}$$

(56)

$${\mathit{\lambda }}_2=\left[\begin{array}{ccc}{\lambda }_2^{11}& {\lambda }_2^{22}& {\lambda }_2^{12}\end{array}\right]$$

(57)

$${\mathit{\lambda }}_3=\left[\begin{array}{cccccc}{\lambda }_3^{111}& {\lambda }_3^{221}& {\lambda }_3^{121}& {\lambda }_3^{112}& {\lambda }_3^{222}& {\lambda }_3^{122}\end{array}\right]$$

(58)

$${\mathit{F}}_2=-{\displaystyle \sum _{e=1}{\displaystyle {\int }_{Y^e}{\mathit{B}}_2^T{\mathit{D}}^{\epsilon }\mathrm{d}Y}}$$

(59)

$${\mathit{F}}_3={\displaystyle \sum _{e=1}{\displaystyle {\int }_{Y^e}\left(2{\mathit{B}}_1^T{\mathit{D}}_1^{\epsilon }+2{\mathit{B}}_1^T\left[\begin{array}{cc}{\mathit{D}}_2^{\epsilon }{\mathit{B}}_2{{\rm X}}_2^e& {\mathit{D}}_3^{\epsilon }{\mathit{B}}_2{{\rm X}}_2^e\end{array}\right]-2{\mathit{B}}_2^T\left[\begin{array}{cc}{\mathit{D}}_4^{\epsilon }{\mathit{B}}_1{{\rm X}}_2^e& {\mathit{D}}_5^{\epsilon }{\mathit{B}}_1{{\rm X}}_2^e\end{array}\right]\right)\mathrm{d}Y}}$$

(60)

$${\mathit{B}}_1=\left[\begin{array}{c}\frac{\partial \mathit{N}}{\partial x_1}\\ \frac{\partial \mathit{N}}{\partial x_2}\end{array}\right], {\mathit{B}}_2=\left[\begin{array}{c}\frac{{\partial }^2\mathit{N}}{\partial x_1^2}\\ \frac{{\partial }^2\mathit{N}}{\partial x_2^2}\\ \frac{2{\partial }^2\mathit{N}}{\partial x_1\partial x_2}\end{array}\right]$$

(61)

$${\mathit{D}}_1^{\epsilon }=\left[\begin{array}{llllll}{D}_{1111}^{\epsilon }\hfill & D_{1122}^{\epsilon }\hfill & D_{1112}^{\epsilon }\hfill & D_{1211}^{\epsilon }\hfill & D_{1222}^{\epsilon }\hfill & D_{1212}^{\epsilon }\hfill \\ D_{1211}^{\epsilon }\hfill & D_{1222}^{\epsilon }\hfill & D_{1212}^{\epsilon }\hfill & D_{2211}^{\epsilon }\hfill & D_{2222}^{\epsilon }\hfill & D_{2212}^{\epsilon }\hfill \end{array}\right]$$

(62)

$${\mathit{D}}_2^{\epsilon }=\left[\begin{array}{lll}{D}_{1111}^{\epsilon }\hfill & D_{1122}^{\epsilon }\hfill & D_{1112}^{\epsilon }\hfill \\ D_{1211}^{\epsilon }\hfill & D_{1222}^{\epsilon }\hfill & D_{1212}^{\epsilon }\hfill \end{array}\right], {\mathit{D}}_3^{\epsilon }=\left[\begin{array}{lll}{D}_{1211}^{\epsilon }\hfill & D_{1222}^{\epsilon }\hfill & D_{1212}^{\epsilon }\hfill \\ D_{2211}^{\epsilon }\hfill & D_{2222}^{\epsilon }\hfill & D_{2212}^{\epsilon }\hfill \end{array}\right]$$

(63)

$${\mathit{D}}_4^{\epsilon }=\left[\begin{array}{ll}{D}_{1111}^{\epsilon }\hfill & D_{1112}^{\epsilon }\hfill \\ D_{2211}^{\epsilon }\hfill & D_{2212}^{\epsilon }\hfill \\ D_{1211}^{\epsilon }\hfill & D_{1212}^{\epsilon }\hfill \end{array}\right], {\mathit{D}}_5^{\epsilon }=\left[\begin{array}{ll}{D}_{1112}^{\epsilon }\hfill & D_{1122}^{\epsilon }\hfill \\ D_{2212}^{\epsilon }\hfill & D_{2222}^{\epsilon }\hfill \\ D_{1212}^{\epsilon }\hfill & D_{1222}^{\epsilon }\hfill \end{array}\right]$$

(64)

where K denotes the global stiffness matrix of the microscopic unit cell model; $Y^e$ is the domain of the element e in the unit cell model; N represents the shape function row vector; ${\mathit{F}}_2$ and ${\mathit{F}}_3$ stand for the second-order and third-order quasi-load matrices, respectively.

Since the node parameters of the plate elements contain deflection and two rotation angles, the periodic conditions (38), (39) and (43) can be implemented by the master-slave method to make the node parameters $\left({\mathit{\chi }}_n,\partial {\mathit{\chi }}_n/\partial y_1,\partial {\mathit{\chi }}_n/\partial y_2\right)$ for n = 2, 3 on slave boundaries be identical to those on the corresponding master boundaries. Then, the nodal degrees of freedom on the slave boundaries are eliminated by applying the multifreedom constraints such that

$${\mathit{X}}_2=\mathit{T}{\tilde{\mathit{X}}}_2$$

(65)

$${\mathit{X}}_3=\mathit{T}{\tilde{\mathit{X}}}_3$$

(66)

where T represents the transformation matrix. By inserting Equations (65) and (66) into Equations (52) and (53), one has

$$\left[\begin{array}{cc}{\mathit{T}}^{\mathrm{T}}\mathit{K}\mathit{T}& {\mathit{T}}^{\mathrm{T}}{\mathit{C}}^{\mathrm{T}}\\ \mathit{C}\mathit{T}& \mathbf{0}\end{array}\right]\left[\begin{array}{c}{\tilde{\mathit{X}}}_2\\ {\lambda }_2\end{array}\right]=\left[\begin{array}{c}{\mathit{T}}^{\mathrm{T}}{\mathit{F}}_2\\ \mathbf{0}\end{array}\right]$$

(67)

$$\left[\begin{array}{cc}{\mathit{T}}^{\mathrm{T}}\mathit{K}\mathit{T}& {\mathit{T}}^{\mathrm{T}}{\mathit{C}}^{\mathrm{T}}\\ \mathit{C}\mathit{T}& \mathbf{0}\end{array}\right]\left[\begin{array}{c}{\tilde{\mathit{X}}}_3\\ {\lambda }_3\end{array}\right]=\left[\begin{array}{c}{\mathit{T}}^{\mathrm{T}}{\mathit{F}}_3\\ \mathbf{0}\end{array}\right]$$

(68)

By solving Equations (67) and (68) in conjunction with Equations (65) and (66), ${\mathit{X}}_2$ and ${\mathit{X}}_3$ are finally achieved. With the obtained ${\mathit{X}}_2$ and Equation (37), the numerical result of the homogenized stiffness matrix of a periodic Kirchhoff plate is achieved as follows:

$${\mathit{D}}^{\mathrm{H}}=\frac{1}{\left|Y\right|}{\displaystyle \sum _{e=1}{\displaystyle {\int }_{Y^e}\left({\mathit{D}}^{\epsilon }+{\mathit{D}}^{\epsilon }{\mathit{B}}_2{\mathit{X}}_2^e\right)\mathrm{d}Y}}$$

(69)

Then for a static problem (14) of the periodic plate, one can easily figure out the homogenized solution by using the finite element methods. For clarity, the finite element formulations for this problem are given as follows:

$${\mathit{K}}^{\mathrm{H}}{\mathit{W}}_0={\mathit{F}}^{\mathrm{H}}$$

(70)

where

$${\mathit{K}}^{\mathrm{H}}={\displaystyle \sum _{e=1}{\displaystyle {\int }_{{\Omega }^e}{\mathit{B}}_2^T{\mathit{D}}^{\mathrm{H}}{\mathit{B}}_2\mathrm{d}\mathrm{\Omega }}}$$

(71)

$${\mathit{F}}^{\mathrm{H}}={\displaystyle \sum _{e=1}{\displaystyle {\int }_{{\Omega }^e}{\mathit{N}}^{T}q\mathrm{d}\mathrm{\Omega }}}$$

(72)

$${\mathit{W}}_0={\displaystyle \sum _{e=1}{\mathit{W}}_0^e}$$

(73)

where ${\mathit{K}}^{\mathrm{H}}$ and ${\mathit{F}}^{\mathrm{H}}$ denote the global stiffness matrix and global load vector of the macroscopic plate model, respectively; ${\mathit{W}}_0^e$ represents the element nodal displacement vector; ${\mathsf{\Omega }}^e$ is the domain of the plate element e.

After solving ${\mathit{W}}_0$ from Equation (70), the asymptotic solution perturbed to the third order in the plate element e can be determined in a matrix form as follows:

$$w_e^{\epsilon }=\mathit{N}{\mathit{W}}_0^e+{\epsilon }^2{\left(\mathit{N}{\mathit{X}}_2^e\right)}^{\mathrm{T}}{\mathit{B}}_2{\mathit{W}}_0^e+{\epsilon }^3{\left(\mathit{N}{\mathit{X}}_3^e\right)}^{\mathrm{T}}{\mathit{B}}_3{\mathit{W}}_0^e$$

(74)

where

$${\mathit{B}}_3=\left[\begin{array}{l}{\partial }^3\mathit{N}/\partial x_1^3\hfill \\ {\partial }^3\mathit{N}/\partial x_1\partial x_2^2\hfill \\ 2{\partial }^3\mathit{N}/\partial x_1^2\partial x_2\hfill \\ {\partial }^3\mathit{N}/\partial x_1^2\partial x_2\hfill \\ {\partial }^3\mathit{N}/\partial x_2^3\hfill \\ 2{\partial }^3\mathit{N}/\partial x_1\partial x_2^2\hfill \end{array}\right]$$

(75)

It follows from Equation (74) that the accuracy of the asymptotic solutions depends on the accuracy of the homogenized displacements and their derivatives if the unit cell problems are solved accurately.

Here, a new method is introduced for calculating higher-order derivatives of the homogenized displacements whose order is beyond the order of the shape functions.

Without loss of generality, take an example of the four-node rectangular element with nodal parameters $\left({\tilde{w}}_n, \partial {\tilde{w}}_n/\partial x_2,-\partial {\tilde{w}}_n/\partial x_1\right)$ $\left(n=1~4\right)$, denoted as ACM12. The form of ${\mathit{W}}_0^e$ containing twelve nodal parameters is as follows:

$${\left({\mathit{W}}_0^e\right)}^{\mathrm{T}}=\left[{\tilde{w}}_{01}^e, \frac{\partial {\tilde{w}}_{01}^e}{\partial x_2}, -\frac{\partial {\tilde{w}}_{01}^e}{\partial x_1}, \cdots , {\tilde{w}}_{04}^e, \frac{\partial {\tilde{w}}_{04}^e}{\partial x_2}, -\frac{\partial {\tilde{w}}_{04}^e}{\partial x_1}\right]$$

(76)

According to the property of the shape function N, one has

$${\tilde{w}}_{0n}^e=\mathit{N}({\mathit{x}}_n){\mathit{W}}_0^e$$

(77)

$$\frac{\partial {\tilde{w}}_{0n}^e}{\partial x_2}=\frac{\partial \mathit{N}({\mathit{x}}_n)}{\partial x_2}{\mathit{W}}_0^e$$

(78)

$$-\frac{\partial {\tilde{w}}_{0n}^e}{\partial x_1}=-\frac{\partial \mathit{N}({\mathit{x}}_n)}{\partial x_1}{\mathit{W}}_0^e$$

(79)

where $\mathit{N}({\mathit{x}}_n)$ represents ${\left.\mathit{N}\right|}_{x=x_n}$, and the same interpretations are also true for other similar expressions. Substituting Equations (77)–(79) into Equation (76) results in:

$$\begin{array}{c}{\left({\mathit{W}}_0^e\right)}^{\mathrm{T}}=\left[\mathit{N}\left({\mathit{x}}_1\right){\mathit{W}}_0^e, \frac{\partial \mathit{N}\left({\mathit{x}}_1\right)}{\partial x_2}{\mathit{W}}_0^e, -\frac{\partial \mathit{N}\left({\mathit{x}}_1\right)}{\partial x_1}{\mathit{W}}_0^e, \cdots ,\right.\\ \left.\mathit{N}\left({\mathit{x}}_4\right){\mathit{W}}_0^e, \frac{\partial \mathit{N}\left({\mathit{x}}_4\right)}{\partial x_2}{\mathit{W}}_0^e, -\frac{\partial \mathit{N}\left({\mathit{x}}_4\right)}{\partial x_1}{\mathit{W}}_0^e\right]\end{array}$$

(80)

Then, the first-order derivative of ${\mathit{w}}_0^e$ can be computed by:

$$\frac{\partial w_0^e}{\partial x_i}=\frac{\partial \mathit{N}}{\partial x_i}{\mathit{W}}_0^e$$

(81)

or

$$\frac{\partial w_0^e}{\partial x_i}=\mathit{N}\frac{\partial {\mathit{W}}_0^e}{\partial x_i}$$

(82)

where

$$\begin{array}{c}{\left(\frac{\partial {\mathit{W}}_0^e}{\partial x_i}\right)}^{\mathrm{T}}=\left[\frac{\partial \mathit{N}\left({\mathit{x}}_1\right)}{\partial x_i}{\mathit{W}}_0^e, \frac{{\partial }^2\mathit{N}\left({\mathit{x}}_1\right)}{\partial x_i\partial x_2}{\mathit{W}}_0^e, -\frac{{\partial }^2\mathit{N}\left({\mathit{x}}_1\right)}{\partial x_i\partial x_1}{\mathit{W}}_0^e, \cdots ,\right.\\ \left.\frac{\partial \mathit{N}\left({\mathit{x}}_4\right)}{\partial x_i}{\mathit{W}}_0^e, \frac{{\partial }^2\mathit{N}\left({\mathit{x}}_4\right)}{\partial x_i\partial x_2}{\mathit{W}}_0^e, -\frac{{\partial }^2\mathit{N}\left({\mathit{x}}_4\right)}{\partial x_i\partial x_1}{\mathit{W}}_0^e\right]\end{array}$$

(83)

Note that above $\partial {\mathit{W}}_0^e/\partial x_i$ is not the direct derivative of the node displacement vector, but the column vector of the nodal displacement partial derivative concerning $x_i$. The same interpretations are also true for the following higher-order terms in Equations (84)–(86).

Replacing ${\mathit{W}}_0^e$ in Equation (81) by $\frac{{\partial }^{r-1}{\mathit{W}}_0^e}{\partial x_j^{r-1}}$ yields:

$$\begin{array}{c}{\left(\frac{\partial {\mathit{W}}_0^e}{\partial x_i}\right)}^{\mathrm{T}}=\left[\frac{\partial \mathit{N}\left({\mathit{x}}_1\right)}{\partial x_i}\frac{{\partial }^{r-1}{\mathit{W}}_0^e}{\partial x_j^{r-1}}, \frac{{\partial }^2\mathit{N}\left({\mathit{x}}_1\right)}{\partial x_i\partial x_2}\frac{{\partial }^{r-1}{\mathit{W}}_0^e}{\partial x_j^{r-1}}, -\frac{{\partial }^2\mathit{N}\left({\mathit{x}}_1\right)}{\partial x_i\partial x_1}\frac{{\partial }^{r-1}{\mathit{W}}_0^e}{\partial x_j^{r-1}}, \cdots ,\right.\\ \left.\frac{\partial \mathit{N}\left({\mathit{x}}_4\right)}{\partial x_i}\frac{{\partial }^{r-1}{\mathit{W}}_0^e}{\partial x_j^{r-1}}, \frac{{\partial }^2\mathit{N}\left({\mathit{x}}_4\right)}{\partial x_i\partial x_2}\frac{{\partial }^{r-1}{\mathit{W}}_0^e}{\partial x_j^{r-1}}, -\frac{{\partial }^2\mathit{N}\left({\mathit{x}}_4\right)}{\partial x_i\partial x_1}\frac{{\partial }^{r-1}{\mathit{W}}_0^e}{\partial x_j^{r-1}}\right]\end{array}$$

(84)

Here, if $r=0$, we define

$$\frac{{\partial }^r{\mathit{W}}_0^e}{\partial x_i^r}={\mathit{W}}_0^e$$

(85)

By using Equations (82)–(85), one can calculate any high-order derivatives of ${\mathit{w}}_0^e$ by:

$$\frac{{\partial }^r{w}_0^e}{\partial x_i^r}=\frac{\partial \mathit{N}}{\partial x_i^{}}\frac{{\partial }^{r-1}{\mathit{W}}_0^e}{\partial x_i^{r-1}}$$

(86)

4. Physical Interpretation of AHM

For a better understanding of the present work, this section aims to interpret the AHM for plates from a physical point of view. In the AHM, the actual (or two-scale) displacements are the superposition of homogenized displacement ${\mathit{w}}_0$ and the perturbed displacements ${\mathit{w}}_n(\mathit{x}, \mathit{y})$, which are the products of ${\mathit{\chi }}_n$ and the n-th derivatives of ${\mathit{w}}_0$. Therefore, the physical interpretation [17] of the AHM is essentially equivalent to that of the influence function ${\mathit{\chi }}_n$, or more specifically, the quasi-load ${\mathit{F}}_n$ for solving ${\mathit{\chi }}_n$.

It can be seen from the unit cell problems (32) and (34) as well as their finite element formulations (59) and (60) that ${\mathit{F}}_n$ (n = 2, 3) completely depend on the material constants and also the microstructure of the unit cell. In other words, ${\mathit{F}}_n$, independent of external forces, are caused by the inhomogeneity of the unit cells. Hence, ${\mathit{F}}_n$ is called quasi-load in the present work, and its values are zero for homogenized structures. Note that ${\mathit{F}}_n$ is self-balanced in a unit cell, implying the zero mean value of ${\mathit{F}}_n$ in Y.

For clarity, a unit cell with a single inclusion, as shown in Figure 5, is used to demonstrate the physical interpretation of the AHM. The unit cell is $1 \mathrm{cm}$ in height, and its in-plane size is $a_1\times a_2=3 \mathrm{cm}\times 3 \mathrm{cm}$ with a square inclusion at its center. Both the matrix and the inclusion are isotropic, with the same Poisson’s ratio of 0.3, and Young’s moduli $E_{\mathrm{I}}=10 \mathrm{GPa}$ for inclusion and $E_{\mathrm{M}}=100 \mathrm{GPa}$ for the matrix, respectively. The unit cell is modeled by $15\times 15$ four-node cubic Hermite rectangular elements here.

Figure 6 and Figure 7, respectively, present the diagrams of each column of ${\mathit{F}}_2$ and ${\mathit{F}}_3$, where Q₃ denotes the nodal shear forces along the thickness direction, M₂₂ and M₁₁ represent the nodal bending moments around the y₁ and y₂ axes, respectively. Some conclusions can be drawn from Figure 6 and Figure 7, as follows:

(1)

The quasi-shear forces are zero for ${\mathit{F}}_2^{11}$ and ${\mathit{F}}_2^{22}$, and the quasi-bending moments for ${\mathit{F}}_2^{12}$ is equal to zero. Compared with ${\mathit{F}}_2$, ${\mathit{F}}_3$ has nonzero quasi-bending moments and nonzero quasi-shear forces, and their distributions are more complex than those of ${\mathit{F}}_2$.

(2)

${\mathit{F}}_2$ is self-balanced in the unit cell, and has nonzero values only at the interfaces of the matrix and inclusion as well as the boundaries, implying that ${\mathit{F}}_2$ is caused by the discontinuities of materials. Additionally, the simple distributions of ${\mathit{F}}_2$ shows that ${\mathit{F}}_2$ is the most fundamental quasi-load reflecting the inhomogeneity of unit cells, implying that ${\mathit{w}}_2$ is the primary perturbed term. Since ${\mathit{F}}_3$ at all nodes is not zero, ${\mathit{F}}_3$ captures more microscopic information compared with ${\mathit{F}}_2$. In general, the AHM perturbed to the third order is accurate enough.

(3)

The unit of ${\mathit{F}}_2$ is $\mathrm{N}·\mathrm{m}$, demonstrating that ${\mathit{F}}_2$ behaves as a moment. In fact, the first two columns of ${\mathit{F}}_2$ are the quasi-bending moments caused by the unit bending curvatures around the ${\mathit{y}}_2$ and ${\mathit{y}}_1$ axes, and the third column of ${\mathit{F}}_2$ denotes the quasi-torsional moment attributable to the unit torsional strain. Accordingly, ${\mathit{\chi }}_2$ represent three fundamental deformations caused by ${\mathit{F}}_2$. Since ${\mathit{w}}_2$ is the product of ${\mathit{\chi }}_2$ and the second derivative of ${\mathit{w}}_0$, see Equation (31), thus, the second derivative of ${\mathit{w}}_0$ acts as the modal coordinates in the superposition method.

(4)

The unit of ${\mathit{F}}_3$ is $\mathrm{N}·{\mathrm{m}}^2$. Since the six columns of ${\mathit{F}}_3$ are independent, thus, ${\mathit{\chi }}_3$ denotes six independent microscopic deformations accordingly, and the third derivative of ${\mathit{w}}_0$ acts as the modal coordinates.

5. Numerical Examples

This section uses several dynamic and static problems of periodic composite plates with different boundary conditions to illustrate the validity and accuracy of the present method.

5.1. Free Vibrations

Two free vibration problems are employed to show the accuracy of the present homogenized elastic constants and the applicability of the Kirchhoff plate theory in the analysis of periodic composite thin plates. For comparison, the 3D FEM results are taken as reference solutions, which are obtained with ABAQUS, and the twenty-node quadratic hexahedral element (C3D20R) is employed. The ACM12 element is used for 2D plates and also unit cell problems.

• Example 1: Free Vibration of a CCCC Single-Layer Periodic Plate

Consider a fully clamped single-layer plate with a size of $L_1\times L_2\times h=30 \mathrm{cm}\times 30 \mathrm{cm}\times 1 \mathrm{cm}$, as shown in Figure 8. The plate with in-plane periodicity consists of two isotropic materials, where Young’s moduli and Poisson’s ratio of the matrix and soft inclusion are the same as those in Section 4, and the mass densities are ${\rho }_{\mathrm{I}}=1142 \mathrm{kg}/{\mathrm{m}}^3$ and ${\rho }_{\mathrm{M}}=2774 \mathrm{kg}/{\mathrm{m}}^3$, respectively.

According to Equation (9), one can obtain the stiffness matrix of the inclusion and the matrix as follows.

$${\mathit{D}}_{\mathrm{I}}^{\epsilon }=\left[\begin{array}{ccc}9.1575& 2.7473& 0\\ 2.7473& 9.1575& 0\\ 0& 0& 3.2051\end{array}\right]\times {10}^2\mathrm{N}\cdot \mathrm{m}$$

(87)

$${\mathit{D}}_{\mathrm{M}}^{\epsilon }=\left[\begin{array}{ccc}9.1575& 2.7473& 0\\ 2.7473& 9.1575& 0\\ 0& 0& 3.2051\end{array}\right]\times {10}^3\mathrm{N}\cdot \mathrm{m}$$

(88)

In this case, the 2D unit cell model is divided into $30\times 30$ ACM12 elements. Figure 9 and Figure 10 present the deformations of the ${\mathit{\chi }}_2$ and ${\mathit{\chi }}_3$ of the periodic plate, respectively. It can be seen that ${\mathit{\chi }}_3$ captures more microscopic properties compared to ${\mathit{\chi }}_2$, refer to the interpretations in Section 4.

The homogenized stiffness matrix ${\mathit{D}}^{\mathrm{H}}$ can be obtained by Equation (69) as:

$${\mathit{D}}^{\mathrm{H}}=\left[\begin{array}{ccc}7.4706& 1.8731& 0\\ 1.8731& 7.4706& 0\\ 0& 0& 2.8352\end{array}\right]\times {10}^3\mathrm{N}\cdot \mathrm{m}$$

(89)

The homogenized density ${\rho }^{\mathrm{H}}$ is determined by the volume average of both materials over the 3D unit cell domain V as:

$${\rho }^{\mathrm{H}}=\frac{1}{\left|V\right|}{\displaystyle {\int }_V{\rho }^{\epsilon }dV}=2.5927\times {10}^3\mathrm{kg}/{\mathrm{m}}^3$$

(90)

Table 1. First five-order natural frequencies of the CCCC single-layer periodic plate.

Order	Frequency/Hz (Relative Error)
	3D FEM	2D FEM	HOM
1	$1.0610\times {10}^3$	$1.0861\times {10}^3 \left(2.36\%\right)$	$1.0787\times {10}^3 \left(1.67\%\right)$
2	$2.1466\times {10}^3$	$2.2181\times {10}^3 \left(3.33\%\right)$	$2.1981\times {10}^3 \left(2.40\%\right)$
3	$2.1466\times {10}^3$	$2.2181\times {10}^3 \left(3.33\%\right)$	$2.1981\times {10}^3 \left(2.40\%\right)$
4	$3.1377\times {10}^3$	$3.2773\times {10}^3 \left(4.45\%\right)$	$3.2296\times {10}^3 \left(2.93\%\right)$
5	$3.8046\times {10}^3$	$3.9821\times {10}^3 \left(4.67\%\right)$	$3.9392\times {10}^3 \left(3.54\%\right)$

and

Table 2. First five-order mode shapes of the CCCC single-layer periodic plate.

Order	3D FEM	HOM
1
2
3
4
5

show the first five-order natural frequencies and mode shapes, respectively, where relative error = (other method—3D FEM)/3D FEM × 100%. Hereinafter, ‘HOM’ represents the homogenized results obtained by the AHM, ‘3D FEM’ stands for the reference results via ABAQUS, and ‘2D FEM’ denotes the FEM’s results based on the Kirchhoff plate theory. The 3D periodic plate is modeled by 24,300 C3D20R elements, while the 2D periodic and homogenized plate models are discretized by 22,500 and 400 ACM12 elements, respectively. It can be seen that the homogenized eigen solutions based on the Kirchhoff plate theory agree well with the 2D and 3D FEM’s results.

• Example 2: Free vibration of an SSSS four-layer periodic plate

Consider a four-layer periodic plate with a size $L_1\times L_2\times h=30 \mathrm{cm}\times 30 \mathrm{cm}\times 1 \mathrm{cm}$. The plate consists of $10\times 10$ unit cells with a size $a_1\times a_2=3 \mathrm{cm}\times 3 \mathrm{cm}$, as shown in Figure 1. Material properties are listed as follows:

• Inclusion: Young’s modulus $E_{\mathrm{I}}=10 \mathrm{GPa}$; Poisson’s ratio $v_{\mathrm{I}}=0.3$; Density ${\rho }_{\mathrm{I}}=1958 \mathrm{kg}/{\mathrm{m}}^3$.
• Matrix: Young’s modulus $E_{\mathrm{M}}=100 \mathrm{GPa}$; Poisson’s ratio $v_{\mathrm{M}}=0.3$; Density ${\rho }_{\mathrm{M}}=2366 \mathrm{kg}/{\mathrm{m}}^3$.

In accordance with Equation (9), one can obtain the stiffness matrix of materials 1 and 2 as follows:

$${\mathit{D}}_{\mathrm{I}}^{\epsilon }=\left[\begin{array}{ccc}1.9460& 0.5838& 0\\ 0.5838& 1.9460& 0\\ 0& 0& 0.6811\end{array}\right]\times {10}^3 \mathrm{N}\cdot \mathrm{m}$$

(91)

$${\mathit{D}}_{\mathrm{M}}^{\epsilon }=\left[\begin{array}{ccc}8.1273& 2.4832& 0\\ 2.4832& 8.1273& 0\\ 0& 0& 2.8446\end{array}\right]\times {10}^3 \mathrm{N}\cdot \mathrm{m}$$

(92)

In this case, the 2D unit cell is discretized into $30\times 30$ ACM12 elements for achieving ${\mathit{\chi }}_2$, ${\mathit{\chi }}_3,$ and ${\mathit{D}}^{\mathrm{H}}$, of which:

$${\mathit{D}}^{\mathrm{H}}=\left[\begin{array}{ccc}4.5384& 1.0667& 0\\ 1.0667& 4.5384& 0\\ 0& 0& 1.8528\end{array}\right]\times {10}^3 \mathrm{N}\cdot \mathrm{m}$$

(93)

and ${\rho }^{\mathrm{H}}$ is calculated as:

$${\rho }^{\mathrm{H}}=\frac{1}{\left|V\right|}{\displaystyle {\int }_V{\rho }^{\epsilon }dV}=2162 \mathrm{kg}/{\mathrm{m}}^3$$

(94)

To further validate the effective properties obtained by the AHM, the free vibration problem of the fully simply supported (SSSS) plate is investigated.

Table 3. First five-order natural frequencies of the SSSS four-layer periodic plate.

Order	Frequency/Hz (Relative Error)
	3D FEM	2D FEM	HOM
1	$4.9631\times {10}^2$	$5.0970\times {10}^2 \left(2.70\%\right)$	$5.1159\times {10}^2 \left(3.08\%\right)$
2	$1.2399\times {10}^3$	$1.2682\times {10}^3 \left(2.28\%\right)$	$1.2722\times {10}^3 \left(2.61\%\right)$
3	$1.2399\times {10}^3$	$1.2682\times {10}^3 \left(2.28\%\right)$	$1.2722\times {10}^3 \left(2.61\%\right)$
4	$1.9512\times {10}^3$	$2.0382\times {10}^3 \left(4.46\%\right)$	$2.0389\times {10}^3 \left(4.49\%\right)$
5	$2.4534\times {10}^3$	$2.5262\times {10}^3 \left(2.97\%\right)$	$2.5349\times {10}^3 \left(3.32\%\right)$

and

Table 4. First five order mode shapes of the SSSS four-layer periodic plate.

Order	3D FEM	HOM
1
2
3
4
5

provide the first five-order natural frequencies and the corresponding mode shapes, respectively, where the 3D model consists of 352,800 C3D20R elements, and the 2D periodic and homogenized plates are discretized by 22,500 and 400 ACM12 elements, respectively.

It can be observed from Table 3 that the frequencies of the 2D models (2D FEM and HOM) match well with those of the 3D models. Additionally, the first five-order mode shapes of the 2D and 3D models shown in Table 4 agree well, implying that the Kirchhoff plate theory is physically suitable for describing the transverse deformation of periodic composite thin plates. Hence, in the following static problems, the 2D FEM displacements obtained with ACM12 elements are regarded as the reference to validate the accuracy of the AHM.

5.2. Static Problems

Two static problems are investigated to further show the accuracy of the homogenized elastic constants and the importance of higher-order perturbed terms in the two-scale asymptotic solutions. In all comparisons, ‘FEM’ represents the referenced results achieved by 2D FEM based on the Kirchhoff plate theory, and ‘AHM2’ and ‘AHM3’ denote the AHM’s results perturbed up to the second and third orders, respectively. For easier comparison, nodes by nodes, ${\mathrm{w}}_0$ and its derivatives are calculated with the same fine grids as those for obtaining referenced solutions or FEM’s solutions with fine meshes, but ${\mathrm{w}}_0$ and the derivatives with the same accuracy can be obtained with coarse grids. Additionally, in order to evaluate the accuracy of the AHM’s results for different orders, a dimensionless residual sum of squares (DRSS) is defined as:

$$\mathrm{DRSS}=\frac{\sqrt{{\displaystyle {\int }_{\Omega }{\left(\mathrm{AHM}-\mathrm{Ref}\right)}^2\mathrm{d}\Omega }}}{{\max }_{\Omega }\left(\mathrm{Ref}\right)}$$

(95)

• Example 3: Static problem of a CFCF periodic layered plate

This example considers the static problem of the periodic layered plate as shown in Figure 4, whose material properties remain constant in the x₂ direction. This layered plate has two clamped opposite sides ($x_1=0 \mathrm{cm}$ and $x_1=30 \mathrm{cm}$) and the other two sides are free, denoted by CFCF here. In addition, this plate shares the same size and material parameters as those in Example 1.

The 2D unit cell of $3 \mathrm{cm}\times 3 \mathrm{cm}$ is also divided into $30\times 30$ ACM12 elements. Figure 11 and Figure 12 present ${\chi }_2$ and ${\chi }_3$ of the periodic layered plate, where ${\chi }_2^{12}$, ${\chi }_3^{121}$, ${\chi }_3^{112},$ and ${\chi }_3^{222}$ are zero. It follows from Figure 11 and Figure 12 that ${\chi }_2$ and ${\chi }_3$ are independent of $y_2$, which is consistent with Equation (45).

After achieving ${\mathit{\chi }}_2$, ${\mathit{D}}^{\mathrm{H}}$ is obtained with Equation (37), as:

$${\mathit{D}}^{\mathrm{H}}=\left[\begin{array}{ccc}2.2894& 0.6868& 0\\ 0.6868& 6.0394& 0\\ 0& 0& 2.2436\end{array}\right]\times {10}^3 \mathrm{N}\cdot \mathrm{m}$$

(96)

Note that ${\mathit{D}}^{\mathrm{H}}$ in Equation (96) is completely the same as the analytical result obtained with Equation (49). Here, the CFCF layered plate is subjected to a uniform transversely distributed load $q={10}^5 \mathrm{N}/{\mathrm{m}}^2$. Figure 13 compares the deflection w and its derivative $\partial w/\partial x_1$ of Line A, and $\partial w/\partial x_2$ is not considered because of the symmetry of the load distribution and structural configuration. Hereinafter, the error in comparisons is defined as follows:

$$\mathrm{Error}=\frac{\left(\mathrm{AHM}-\mathrm{Ref}\right)}{{\max }_{Line}\left(\mathrm{Ref}\right)}$$

(97)

It can be seen from Figure 13 that w and $\partial w/\partial x_1$ of the AHM match well with the referenced results. Additionally, it can also be observed from Figure 13 that the results of both AHM2 and AHM3 do not satisfy the displacement boundary conditions, since the influence functions are solved with periodic boundary conditions rather than the plate boundary conditions. In order to evaluate the global accuracy of the AHM’s results,

Table 5. DRSS for the CFCF periodic layered plate.

Perturbed Terms	$\mathit{w}$	$\partial \mathit{w}/\partial {\mathit{x}}_1$
HOM	2.6404	10.0475
AHM2	2.5986	2.7112
AHM3	2.6003	1.4978

compares the DRSS of w and $\partial w/\partial x_1$. It can be concluded from Table 5 that there are no significant differences between w achieved by the HOM, AHM2, and AHM3, but apparently, the higher the perturbation order, the higher the accuracy of $\partial w/\partial x_1.$ This is because higher-order perturbation terms can capture more information about microstructures, thus, AHM2 and AHM3 can achieve more accurate results than HOM.

• Example 4: Static problem of an SFSF four-layer periodic plate

This example shows the static behavior of the SFSF four-layer periodic plate with two opposite edges ($x_2=0$ and $x_2=L_2$) simply supported and the other two edges free, see Figure 1, of which the free vibration has been discussed in Example 2. Here, the plate is subjected to a distributed load $q={10}^7\sin \left(3\mathsf{\pi }{x}_1/L_1\right)\sin \left(4\mathsf{\pi }{x}_2/L_2\right)\mathrm{N}/{\mathrm{m}}^2$.

Figure 14 and Figure 15 compare the displacements along Line A and Line B, respectively.

Table 6. DRSS for the SFSF four-layer periodic plate.

Perturbed Terms	$\mathit{w}$	$\partial \mathit{w}/\partial {\mathit{x}}_2$	$\partial \mathit{w}/\partial {\mathit{x}}_1$
HOM	0.6897	4.1701	2.5050
AHM2	0.5457	2.0187	1.8187
AHM3	0.5603	1.8169	1.6351

shows the DRSS results for w, $\partial w/\partial x_1$ and $\partial w/\partial x_2$ up to different orders. The same conclusions as those in Example 3 can be achieved. It follows from Figure 14 and Figure 15 and Table 6 that the results of AHM2 and AHM3 match better than the homogenized solutions, especially for the rotation angles ($\partial w/\partial x_1$ and $\partial w/\partial x_2$). In addition, it is worth noting that as the power of the external distributed load function increases, the deformations of the plate become more complex, and thus, the second-order or even higher-order perturbed displacements are more necessary for ensuring accuracy.

6. Conclusions

This paper proposes a two-scale method by combining the Kirchhoff plate theory with the two-scale asymptotic homogenization method to deal with the static and dynamic problems of 3D periodic thin plates. In this work, the solutions of the fourth-order elliptic PDE with periodically oscillating coefficients were given in an asymptotic expansion form, where perturbed terms were the multiplications of influence functions and the derivatives of homogenized displacements. To determine the influence functions from the unit cell problems, periodic boundary and normalization constraint conditions were given and elaborated physically.

It was found that the first-order perturbed term of the asymptotic expansion solution of the fourth-order PDE with periodically oscillating coefficients should be zero. The reason for this phenomenon is that the first-order perturbed terms cannot reflect the microdeformations. In addition, with the physical interpretation of the AHM, it was shown that the second-order influence functions are the fundamental terms for capturing the microscopic information, since the second-order quasi-loads for solving the second-order influence functions are simple line loads, with nonzero values only along the interfaces of matrices and inclusions as well as boundaries.

Finally, the free vibrations and static problems of several periodic composite plates with different boundary conditions were investigated to validate the effectiveness of the proposed method, showing that the present AHM’s solutions are meaningful and physically acceptable. Additionally, it has been shown that homogenized displacements play a significant part in the prediction of microscale solutions and that second-order or even higher-order perturbed displacements are necessary for achieving accurate rotation angles for periodic plates. This work lays the foundation for the study of moderately thick periodic plates.