Structural Analysis of Composite Conical Convex-Concave Plate (CCCP) Using VAM-Based Equivalent Model

Abstract

Compared with the ordinary foundation plate, the composite conical convex-concave plate (CCCP) has obvious anisotropic characteristics, and there is less research on the relationship between its mechanical properties and structural parameters. In this article, a numerical model for the equivalent stiffness of a typical unit cell with conical convex is established by using the variational asymptotic method. Then, the 3D finite element model (3D-FEM) of CCCP is transformed into 2D equivalent plate model (2D-EPM) with the effective plate properties obtained from the constitutive analysis of unit cell. The accuracy of 2D-EPM is verified by comparing with the displacement, natural frequencies, and buckling results from 3D-FEM under different boundary conditions. Then, the influence of geometric parameters and layup configurations on the effective performances of CCCP are investigated. Finally, the buckling loads and natural frequencies of bidirectional CCCP are compared with those of CCCP by using the present model. The present model is particularly useful in the early design stage of CCCP where many design trade-offs need to be made over a vast design space in terms of material selection, ply angles, and geometric parameters.

1. Introduction

With the development of modern science and technology, thin plate theory cannot meet the requirement of industrial production. The main disadvantage is that the increase in stiffness can only be achieved by changing the thickness of thin plate, but the thin plate theory is no longer applicable when the thickness of thin plate reaches a certain value, and the shear effect needs to be considered. In addition, the increase in stiffness can also be achieved by bonding laminated plates and sandwich plates, which is not suitable for industrial production and cannot meet engineering requirements [1,2]. Moreover, subsequent studies have found that the production process of laminated plates and sandwich plates is relatively complex and prone to degumming or dislocation, which is not conducive to the use in industrial production [3,4,5].

As a new type of lightweight structure, the convex-concave plate, made by stamping and rolling can avoid the shortcomings of the degumming, dislocation, and complicated process, resulting in high industrial application value and reliability. Although its thickness is very small, it has light weight, high strength, and high stiffness. This is because the convexities of the plate changes the plate structure after stamping or rolling, and the stiffness is greatly enhanced. Furthermore, the convexities can absorb a lot of energy in the process of impact vibration, and the loss of the plate can be greatly reduced [6]. This kind of plate also has high heat dissipation efficiency and excellent sound insulation and heat insulation performance, which cannot be achieved by the sandwich plate and laminated plate. In fact, it has been used as the plate heat exchanger as shown in Figure 1.

Due to the existence of the convexities, the energy absorption capacity is obviously improved and the service life is greatly enhanced [7,8,9]. Haldar et al. [10] focused on studying the influence of the thickness of convexities on the energy absorption and obtained the optimal energy absorption properties of the convex-concave plate by comparing the experimental and finite element simulation results. Sashikumar [11] tested the energy absorption capacity of an aluminum egg-box plate in commercial vehicles and proposed that the natural constraints within the geometric structure of egg-box plates have a positive impact on their energy absorption capacity. Akisanya and Fleck [12] discussed the dependence of strength and energy absorption of egg-box material on the geometric shape by using the finite element prediction method and determined the variation of the stiffness and strength of egg-box material with relative density. Ma et al. [13] concluded that the periodic origami structure has better energy absorption capacity by comparing the dynamic analysis and experimental and numerical results.

In addition to the traditional single-layer convex-concave plate, the laminated convex-concave plate has also been widely investigated. Chang et al. [14,15,16] adopted the experimental method to study the energy absorbing ability and deformation process of two kinds of sinusoidal convex-concave plate reinforced by fiber materials during quasi-static compression. Zhang et al. [17] tested the performance of egg-box panel stuffed whipple shield from impact area, cell size and the axial offset etc. At present, experiments and finite element simulation are mostly used to investigate the energy absorption capacity, failure, and heat transfer performance of concave-convex plates.

The laminated composite plate by nature has two dimensions larger than the thickness dimension by an order of magnitude. In recent years, the variational asymptotic method (VAM) has been used to strictly split the plate/shell problem into the through-the-thickness linear analysis (1-D analysis) and 2D nonlinear plate analysis by using the small parameter of width-to-thickness ratio. This method combines advantages of both asymptotic and variational methods, and considers all possible deformations during problem formulations while avoiding any kinematic assumptions. Atilgan and Hodges [18,19] developed VAM for laminated plates in which each lamina exhibits monoclinic symmetry about its own midplane. Sutyrin and Hodges [20] extended VAM to laminated plate, the material properties of which varied through the thickness and for which each lamina was orthotropic. Later, Sutyrin [21] developed linear, asymptotically correct theories for inhomogeneous orthotropic plates, for example, laminated plates with orthotropic laminae. Yu et al. [22] developed an accurate stress/strain recovery procedure for laminated plates that can be implemented in standard finite element programs. Kamineni and Burela [23] developed constraint method for laminated composite flat stiffened panel using VAM. Further, Zhong and Yu [24,25,26,27,28] successfully constructed the equivalent model of various composite structures by using this method.

The composite conical convex-concave plate (CCCP) studied in this article is a new type of convex-concave plate. Most of the research at this stage is limited to the manufacturing process, and its equivalent mechanical properties have not been investigated, especially the natural vibration and stability characteristics. Therefore, it is necessary to study the equivalent stiffness of the orthotropic convex-concave plate, which lays a foundation for studying the deformation of convex-concave plates under applied load and provides a theoretical basis for the reasonable design of its geometric parameters.

In this work, a VAM-based equivalent plate model is developed to predict the effective performance of CCCP. The organization of this paper is as follows. The theoretical formulation for the VAM-based 2D equivalent plate model (2D-EPM) is deduced in Section 2. Numerical examples of static displacement, global buckling, and free-vibration analysis of CCCP are used in Section 3 to verify the accuracy and effectiveness of 2D-EPM. Section 4 investigates the influences of layup configurations and geometric parameters on the effective performances of CCCP by using the 2D-EPM. In Section 5, the effective performance of bidirectional CCCP and CCCP is compared. Finally, some conclusions are drawn in Section 6.

2. Theoretical Foundation

2.1. Equivalent Plate Modeling of CCCP

The equivalent model of CCCP can be established from the perspective of energy concept. Due to the periodicity of the CCCP along the two axes of the plane, a square plate with conical convexity is taken as a typical unit cell as shown in Figure 2b. The length, height, and thickness of the unit cell are l, h, and t, respectively. That is, the analysis of the original CCCP is decomposed into a microscopic analysis of the unit cell (providing effective plate properties) and macro analysis of 2D-EPM as shown in Figure 2c.

To facilitate the derivation, two groups of coordinates are introduced: the macro coordinates, $x_i$, describing the original structure, and the micro coordinates, $y_i$, describing the unit cell. The macro coordinate $x_3$ is perpendicular to the in-plane coordinates $x_1$ and $x_2$, and the origin of the macro coordinate is located at the center of the plate as shown in Figure 2a. Because the micro size of the unit cell is much smaller than the macro size of the plate, the micro coordinates $y_i=x_i/\xi$ ($1/\xi$ is a small parameter) are used to describe unit cell. For the equivalent plate model, the function of the original CCCP can be expressed as a function defined along the reference plane $x_1-x_2$ ($x_3$ disappears), and its partial derivative is

$$\frac{\partial f\left(x_{\alpha };y_i\right)}{\partial x_{\alpha }}={\left.\frac{\partial f\left(x_{\alpha };y_i\right)}{\partial x_{\alpha }}\right|}_{y_i=const}+{\left.\frac{1}{\xi }\frac{\partial f\left(x_{\alpha };y_i\right)}{\partial y_i}\right|}_{x_{\alpha }=const}\equiv f_{,\alpha }+\frac{1}{\xi }{f}_{;i}$$

(1)

where i takes values of 1, 2, and 3; $\alpha$ takes values 1 and 2.

To establish the VAM-based equivalent model of CCCP, the 2D plate variables should be used to represent the 3D displacements of the original CCCP as

$$\begin{array}{c}{u}_1\left(x_{\alpha };y_i\right)=\underline{v_1\left(x_1,x_2\right)-\xi y_3{v}_{3,1}\left(x_1,x_2\right)}+\xi w_1\left(x_{\alpha };y_i\right)\hfill \\ u_2\left(x_{\alpha };y_i\right)=\underline{v_2\left(x_1,x_2\right)-\xi y_3{v}_{3,2}\left(x_1,x_2\right)}+\xi w_2\left(x_{\alpha };y_i\right)\hfill \\ u_3\left(x_{\alpha };y_i\right)=\underline{v_3\left(x_1,x_2\right)}+\xi w_3\left(x_{\alpha };y_i\right)\hfill \end{array}$$

(2)

where $u_i$ and $v_i$ represent the displacement of the original three-dimensional plate and two-dimensional equivalent plate, respectively, and $w_i$ denotes the fluctuating function used to describe the deformation that cannot be described by the classical plate theory.

The underlined terms in Equation (2) are the deformations generated by the reference plane, which should meet the following requirements:

$$h{v}_{\alpha }\left(x_{\alpha }\right)=〈u_{\alpha }〉+〈\xi y_3〉v_{3,\alpha },h{v}_3\left(x_{\alpha }\right)=〈u_3〉$$

(3)

where the angle brackets denote the integral over the unit cell.

The definitions in Equation (3) introduce the constraints on the fluctuating functions as

$$〈\xi w_i〉=0$$

(4)

The strain field can be obtained from linear elasticity theory as

$${\mathrm{\Gamma }}_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)$$

(5)

The corresponding 3D linear strain field can be obtained by substituting Equation (2) into Equation (5), and neglecting higher-order terms,

$$\begin{array}{c}{\mathrm{\Gamma }}_{11}={\epsilon }_{11}+\xi y_3{\kappa }_{11}+w_{1,1}\hfill \\ 2{\mathrm{\Gamma }}_{12}=2{\epsilon }_{22}+2\xi y_3{\kappa }_{12}+w_{1,2}+w_{2,1}\hfill \\ {\mathrm{\Gamma }}_{22}={\epsilon }_{22}+\xi y_3{\kappa }_{22}+w_{2,2}\hfill \\ 2{\mathrm{\Gamma }}_{13}=w_{1,3}+w_{3,1}\hfill \\ 2{\mathrm{\Gamma }}_{23}=w_{2,3}+w_{3,2}\hfill \\ {\mathrm{\Gamma }}_{33}=w_{3,3}\hfill \end{array}$$

(6)

where the 2D plate strains ${\epsilon }_{\alpha \beta }$ and ${\kappa }_{\alpha \beta }$ of the 2D-EPM are defined as

$${\epsilon }_{\alpha \beta }\left(x_1,x_2\right)=\frac{1}{2}\left(v_{\alpha ,\beta }+v_{\beta ,\alpha }\right),{\kappa }_{\alpha \beta }\left(x_1,x_2\right)=-v_{3,\alpha \beta }$$

(7)

With the aid of matrix representation, the 3D strain field can be expressed as

$$\begin{array}{c}{\mathbf{\Gamma }}_e={\left[{\mathrm{\Gamma }}_{11}{\mathrm{\Gamma }}_{22}2{\mathrm{\Gamma }}_{12}\right]}^{\mathrm{T}}=\mathbf{\epsilon }+x_3\mathbf{\kappa }+{\partial }_e{\mathit{w}}_{\left|\right|}\hfill \\ 2{\mathbf{\Gamma }}_s={\left[2{\mathrm{\Gamma }}_{13}2{\mathrm{\Gamma }}_{23}\right]}^{\mathrm{T}}={\mathit{w}}_{\left|\right|}+{\partial }_t{w}_3\hfill \\ {\mathbf{\Gamma }}_t={\mathrm{\Gamma }}_{33}=w_{3,3}\hfill \end{array}$$

(8)

where ${\mathbf{\Gamma }}_e,{\mathbf{\Gamma }}_s,{\mathbf{\Gamma }}_t$ are strain components of 3D-FEM, respectively; ${\left(\right)}_{\left|\right|}={\left[{\left(\right)}_1{\left(\right)}_2\right]}^{\mathrm{T}}$, $\mathbf{\epsilon }={\left[{\epsilon }_{11}2{\epsilon }_{12}{\epsilon }_{22}\right]}^{\mathrm{T}}$, $\mathbf{\kappa }={\left[{\kappa }_{11}{\kappa }_{12}+{\kappa }_{21}{\kappa }_{22}\right]}^{\mathrm{T}}$, and

$${\partial }_e=\left[\begin{array}{cc}{\left(\right)}_{,1}& 0\\ {\left(\right)}_{,2}& {\left(\right)}_{,1}\\ 0& {\left(\right)}_{,2}\end{array}\right],{\partial }_t=\left\{\begin{array}{c}{\left(\right)}_{,1}\hfill \\ {\left(\right)}_{,2}\hfill \end{array}\right\}$$

(9)

The geometric structure of the element can be divided into five parts shown in Figure 3 for easy integration, and the strain energy of unit cell can be written as

$$\begin{array}{c}{\displaystyle U={\int }_{-l_3/2}^{l_3/2}{\int }_{-\left(l_3/2+l_1+l_2\right)}^{l_3/2+l_1+l_2}{\int }_0^h{\mathbf{\Gamma }}_C^T{\mathit{D}}_C{\mathbf{\Gamma }}_C\mathrm{d}{y}_1\mathrm{d}{y}_2\mathrm{d}{y}_3}\hfill \\ {\displaystyle +2{\int }_{l_3/2}^{l_3/2+l_2}{\int }_{-\left(l_3/2+l_1+l_2\right)}^{l_3/2+l_1+l_2}{\int }_0^h{\mathbf{\Gamma }}_B^T{\mathit{D}}_B{\mathbf{\Gamma }}_B\mathrm{d}{y}_1\mathrm{d}{y}_2\mathrm{d}{y}_3}\hfill \\ {\displaystyle +2{\int }_{l_3/2+l_2}^{l_3/2+l_2+l_1}{\int }_{-\left(l_3/2+l_1+l_2\right)}^{l_3/2+l_1+l_2}{\int }_{-t}^0{\mathbf{\Gamma }}_A^T{\mathit{D}}_A{\mathbf{\Gamma }}_A\mathrm{d}{y}_1\mathrm{d}{y}_2\mathrm{d}{y}_3}\hfill \end{array}$$

(10)

Equation (10) can be briefly expressed as

$$U=\frac{1}{2}〈{\mathbf{\Gamma }}^{\mathrm{T}}\mathit{D}\mathbf{\Gamma }〉=\frac{1}{2}〈{\left\{\begin{array}{c}{\mathbf{\Gamma }}_e\\ 2{\mathbf{\Gamma }}_s\\ {\mathbf{\Gamma }}_t\end{array}\right\}}^{\mathrm{T}}\left[\begin{array}{ccc}{\mathit{D}}_e& {\mathit{D}}_{es}& {\mathit{D}}_{et}\\ {\mathit{D}}_{es}^{\mathrm{T}}& {\mathit{D}}_s& {\mathit{D}}_{st}\\ {\mathit{D}}_{et}^{\mathrm{T}}& {\mathit{D}}_{st}^{\mathrm{T}}& {\mathit{D}}_t\end{array}\right]\left\{\begin{array}{c}{\mathbf{\Gamma }}_e\\ 2{\mathbf{\Gamma }}_s\\ {\mathbf{\Gamma }}_t\end{array}\right\}〉$$

(11)

where ${\mathit{D}}_e,{\mathit{D}}_{es},{\mathit{D}}_{et},{\mathit{D}}_s,{\mathit{D}}_{st}$, and ${\mathit{D}}_t$ denote the sub-matrices of the three-dimensional $6\times 6$ material matrix.

The virtual work done due to applied loads can be expressed as

$$\delta {\overline{W}}_{3D}=\delta {\overline{W}}_{2D}+\delta {\overline{W}}^{*}$$

(12)

where $\delta {\overline{W}}_{2D}$ and $\delta {\overline{W}}^{*}$ can be defined as

$$\delta {\overline{W}}_{2D}=〈p_i\delta v_i+q_{\alpha }\delta v_{3,\alpha }〉,\delta {\overline{W}}^{*}=〈〈f_i\delta w_i〉+{\tau }_i\delta w_i^{+}+{\beta }_i\delta w_i^{-}〉$$

(13)

where ${(·)}^{+}={\left.(·)\right|}_{x_3=t/2}$ and ${(·)}^{-}={\left.(·)\right|}_{x_3=-t/2}$ denote the items acting on the top and bottom surfaces, respectively; $f_i$ is the body force; ${\tau }_i$ and ${\beta }_i$ denote the traction forces on the top and bottom surface, respectively; and $p_i=〈f_i〉+{\tau }_i+{\beta }_i,q_{\alpha }=h/2\left({\beta }_{\alpha }-{\tau }_{\alpha }\right)-〈x_3{f}_{\alpha }〉$.

The total potential energy of unit cell can be expressed as

$$\delta \mathrm{\Pi }=\delta U-\delta W^{*}=\frac{1}{2}\delta 〈{\mathbf{\Gamma }}^{\mathrm{T}}\mathit{D}\mathbf{\Gamma }〉-〈〈f_i\delta w_i〉+{\tau }_i\delta w_i^{+}+{\beta }_i\delta w_i^{-}〉$$

(14)

2.2. Dimension Reduction of CCCP

To solve the unknown fluctuating function $w_i$ by using VAM, the order of each term in Equation (14) must be evaluated as

$$\begin{array}{c}{\displaystyle {\mathrm{\Gamma }}_{ij}\sim {\epsilon }_{\alpha \beta }\sim h{\kappa }_{\alpha \beta }\sim n,w_i\sim hn,{\mathit{w}}_{\left|\right|,\alpha }\sim w_{3,\alpha }\sim \frac{h}{L}n}\hfill \\ {\displaystyle {\mathit{w}}_{\left|\right|,3}\sim w_{3,3}\sim n,h{f}_{\alpha }\sim {\alpha }_{\alpha }\sim {\beta }_{\alpha }\sim \mu \frac{h}{L}n,h{f}_3\sim {\alpha }_3\sim {\beta }_3\sim \mu {\left(\frac{h}{L}\right)}^2}\hfill \end{array}$$

(15)

where n and $\mu$ are the order of the minimum strain and material properties, respectively, and L is the length of the plate.

2.2.1. Zeroth-Order Approximation

The explicit expression of Equation (14) can be expressed as

$$\begin{array}{c}2\mathrm{\Pi }=〈\underline{\underline{{\left(\mathbf{\epsilon }+x_3\mathbf{\kappa }\right)}^{\mathrm{T}}{\mathit{D}}_e\left(\mathbf{\epsilon }+x_3\mathbf{\kappa }\right)}}+\underline{2{\left(\mathbf{\epsilon }+x_3\mathbf{\kappa }\right)}^{\mathrm{T}}{\mathit{D}}_e{\partial }_e{\mathit{w}}_{\left|\right|,\alpha }}\hfill \\ +\underline{2{\left({\partial }_e{\mathit{w}}_{\left|\right|,\alpha }\right)}^{\mathrm{T}}{\mathit{D}}_e{\partial }_e{\mathit{w}}_{\left|\right|}}+2{\left(\mathbf{\epsilon }+x_3\mathbf{\kappa }\right)}^{\mathrm{T}}{\mathit{D}}_{es}{\mathit{w}}_{\left|\right|,3}+\underline{2{\left(\mathbf{\epsilon }+x_3\mathbf{\kappa }\right)}^{\mathrm{T}}{\mathit{D}}_{es}{\partial }_t{w}_{3,\alpha }}\hfill \\ +\underline{2{\left({\partial }_e{\mathit{w}}_{\left|\right|,\alpha }\right)}^{\mathrm{T}}{\mathit{D}}_{es}\left({\mathit{w}}_{\left|\right|,3}+{\partial }_t{w}_{3,\alpha }\right)}+2{\left(\mathbf{\epsilon }+x_3\mathbf{\kappa }\right)}^{\mathrm{T}}{\mathit{D}}_{et}{w}_{3,3}\hfill \\ +\underline{2{\left({\partial }_e{w}_{\left|\right|,\alpha }\right)}^{\mathrm{T}}{\mathit{D}}_{et}{w}_{3,3}}+{\mathit{w}}_{\left|\right|,3}^{\mathrm{T}}{\mathit{D}}_s{w}_{\left|\right|,3}+\underline{2{\mathit{w}}_{\left|\right|,3}^{\mathrm{T}}{\mathit{D}}_s{\partial }_t{w}_{3,\alpha }}+\underline{2{\left({\partial }_t{w}_{3,\alpha }\right)}^{\mathrm{T}}{\mathit{D}}_s{\partial }_t{w}_{3,\alpha }}\hfill \\ +2{\mathit{w}}_{\left|\right|,3}^{\mathrm{T}}{\mathit{D}}_{st}{w}_{3,3}+\underline{2{\left({\partial }_t{w}_{3,\alpha }\right)}^{\mathrm{T}}{\mathit{D}}_{st}{w}_{3,3}}+{\mathit{D}}_t{w}_{3,3}^2〉\hfill \\ -\underline{2\left(〈f_i^{\mathrm{T}}{w}_i〉+{\tau }_i^{\mathrm{T}}{w}_i^{\mathrm{T}}+{\beta }_i^{\mathrm{T}}{w}_i^{\mathrm{T}}\right)}\hfill \end{array}$$

(16)

The fluctuating function can be constrained by introducing the Lagrangian multiplier ${\lambda }_i$, such as

$$\delta \left(\mathrm{\Pi }+{\lambda }_i〈w_i〉\right)=0$$

(17)

The fluctuating function can be obtained by solving the following zeroth-order approximate variational statement after removing the underline and double underline items from Equation (16),

$$〈\begin{array}{c}\left[{\left(\epsilon +x_3\mathbf{\kappa }\right)}^{\mathrm{T}}{\mathit{D}}_{es}+{\mathit{w}}_{\left|\right|,3}^{\mathrm{T}}{\mathit{D}}_s+w_{3,3}^{\mathrm{T}}{\mathit{D}}_{st}^{\mathrm{T}}\right]\delta {\mathit{w}}_{\left|\right|,3}\hfill \\ +{\lambda }_i\delta w_i+\left[{\left(\mathbf{\epsilon }+x_3\mathbf{\kappa }\right)}^{\mathrm{T}}{\mathit{D}}_{et}+{\mathit{w}}_{\left|\right|,3}^{\mathrm{T}}{\mathit{D}}_{st}+w_{3,3}^{\mathrm{T}}{\mathit{D}}_t\right]\delta w_{3,3}\hfill \end{array}〉=0$$

(18)

This results in the following Euler–Lagrange equations:

$$\begin{array}{c}{\left[{\left(\mathbf{\epsilon }+x_3\mathbf{\kappa }\right)}^{\mathrm{T}}{\mathit{D}}_{es}+{\mathit{w}}_{\left|\right|,3}^{\mathrm{T}}{\mathit{D}}_s+w_{3,3}^{\mathrm{T}}{\mathit{D}}_{st}^{\mathrm{T}}\right]}_{,3}={\mathbf{\lambda }}_{\left|\right|}\hfill \\ {\left[{\left(\mathbf{\epsilon }+x_3\mathbf{\kappa }\right)}^{\mathrm{T}}{\mathit{D}}_{et}+{\mathit{w}}_{\left|\right|,3}^{\mathrm{T}}{\mathit{D}}_{st}+w_{3,3}^{\mathrm{T}}{\mathit{D}}_t\right]}_{,3}={\lambda }_3\hfill \end{array}$$

(19)

where ${\mathbf{\lambda }}_{\left|\right|}={\left[{\lambda }_1{\lambda }_2\right]}^{\mathrm{T}}$ and ${\lambda }_3$ are Lagrange multipliers corresponding to the constraint components of ${\mathit{w}}_{\left|\right|}$ and $w_3$.

The boundary conditions are

$$\begin{array}{c}{\left[{\left(\mathbf{\epsilon }+x_3\mathbf{\kappa }\right)}^{\mathrm{T}}{\mathit{D}}_{es}+{\mathit{w}}_{\left|\right|,3}^{\mathrm{T}}{\mathit{D}}_s+w_{3,3}{\mathit{D}}_{st}^{\mathrm{T}}\right]}^{+/-}=0\hfill \\ {\left[{\left(\mathbf{\epsilon }+x_3\mathbf{\kappa }\right)}^{\mathrm{T}}{\mathit{D}}_{et}+{\mathit{w}}_{\left|\right|,3}^{\mathrm{T}}{\mathit{D}}_{st}+w_{3,3}{\mathit{D}}_t\right]}^{+/-}=0\hfill \end{array}$$

(20)

where the superscript “+/−” denotes the items applied on the top and bottom surface of the plate.

We can solve $w_i$ by substituting Equation (20) back into Equation (19):

$${\mathit{w}}_{\left|\right|}={〈-\left(\mathbf{\epsilon }+x_3\mathbf{\kappa }\right){\overline{\mathit{D}}}_{es}{\mathit{D}}_s^{-1}〉}^{\mathrm{T}},w_3=〈-\left(\mathbf{\epsilon }+x_3\mathbf{\kappa }\right){\overline{\mathit{D}}}_{et}{\mathit{D}}_t^{-1}〉$$

(21)

where

$${\overline{\mathit{D}}}_{es}={\mathit{D}}_{es}-{\overline{\mathit{D}}}_{et}{\mathit{D}}_{st}^T{\overline{\mathit{D}}}_t^{-1},{\overline{\mathit{D}}}_{et}={\mathit{D}}_{et}-{\mathit{D}}_{es}{\mathit{D}}_s^{-1}{\mathit{D}}_{st},{\overline{\mathit{D}}}_t={\mathit{D}}_t-{\mathit{D}}_{sl}^T{\mathit{D}}_s^{-1}{\mathit{D}}_{st}$$

(22)

Substituting Equation (22) into Equation (18), we obtain the first approximation of the strain energy as

$$U_{2D}=\frac{1}{2}〈{\left(\mathbf{\epsilon }+x_3\mathbf{\kappa }\right)}^{\mathrm{T}}{\overline{\mathit{D}}}_e\left(\mathbf{\epsilon }+x_3\mathbf{\kappa }\right)〉=\frac{1}{2}{\left\{\begin{array}{c}\mathbf{\epsilon }\hfill \\ \mathbf{\kappa }\hfill \end{array}\right\}}^{\mathrm{T}}\left[\begin{array}{cc}\mathit{A}& \mathit{B}\\ {\mathit{B}}^T& \mathit{D}\end{array}\right]\left\{\begin{array}{c}\mathbf{\epsilon }\hfill \\ \mathbf{\kappa }\hfill \end{array}\right\}$$

(23)

where

$$\begin{array}{c}\mathit{A}=〈{\overline{\mathit{D}}}_e〉,\mathit{B}=〈x_3{\overline{\mathit{D}}}_e〉,\mathit{D}=〈x_3^2{\overline{\mathit{D}}}_e〉,\hfill \\ {\overline{\mathit{D}}}_e={\mathit{D}}_e-{\overline{\mathit{D}}}_{es}{\mathit{D}}_s^{-1}{\mathit{D}}_{es}^{\mathrm{T}}-{\overline{\mathit{D}}}_{et}{\mathit{D}}_{et}^{\mathrm{T}}/{\overline{\mathit{D}}}_t\hfill \end{array}$$

(24)

For the zeroth-order approximate, the 3D strain field can be recovered by using Equation (6) as

$${\mathbf{\Gamma }}_e^0=\mathbf{\epsilon }+x_3\mathbf{\kappa },2{\mathbf{\Gamma }}_s^0=-{\mathit{w}}_{\left|\right|,3},{\mathbf{\Gamma }}_t^0=w_{3,3}$$

(25)

The local stress field can be recovered as

$$\begin{array}{c}{\mathbf{\sigma }}_e^0={\left[{\sigma }_{11}^0{\sigma }_{12}^0{\sigma }_{22}^0\right]}^T=\overline{{\mathit{D}}_e}\left(\mathbf{\epsilon }+x_3\mathbf{\kappa }\right)\hfill \\ {\mathbf{\sigma }}_s^0={\left[{\sigma }_{13}^0{\sigma }_{23}^0\right]}^T=0\hfill \\ {\mathbf{\sigma }}_t^0={\sigma }_{33}^0=0\hfill \end{array}$$

(26)

The resultant stress of the plate can be defined as

$$\begin{array}{c}{\displaystyle \frac{\partial U_{2D}}{\partial {\epsilon }_{11}}=N_{11},\frac{\partial U_{2D}}{\partial 2{\epsilon }_{12}}=N_{12},\frac{\partial U_{2D}}{\partial {\epsilon }_{22}}=N_{22}}\hfill \\ {\displaystyle \frac{\partial U_{2D}}{\partial {\kappa }_{11}}=M_{11},\frac{\partial U_{2D}}{\partial 2{\kappa }_{12}}=M_{12},\frac{\partial U_{2D}}{\partial {\kappa }_{22}}=M_{22}}\hfill \end{array}$$

(27)

The constitutive relation of CCCP can be obtained by connecting the internal stress, strain, and curvature of the plate as

$$\left\{\begin{array}{c}{N}_{11}\hfill \\ N_{22}\hfill \\ N_{12}\hfill \\ M_{11}\hfill \\ M_{22}\hfill \\ M_{12}\hfill \end{array}\right\}=\left[\begin{array}{cccccc}{A}_{11}& A_{12}& A_{16}& B_{11}& B_{12}& B_{16}\\ A_{12}& A_{22}& A_{26}& B_{12}& B_{22}& B_{26}\\ A_{16}& A_{26}& A_{66}& B_{16}& B_{26}& B_{66}\\ B_{11}& B_{12}& B_{16}& D_{11}& D_{12}& D_{16}\\ B_{12}& B_{22}& B_{26}& D_{12}& D_{22}& D_{26}\\ B_{16}& B_{26}& B_{66}& D_{16}& D_{26}& D_{66}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_{11}\\ {\epsilon }_{22}\\ 2{\epsilon }_{12}\\ {\kappa }_{11}\\ {\kappa }_{22}\\ 2{\kappa }_{12}\end{array}\right\}$$

(28)

2.2.2. First-Order Approximation

It can be seen from Equation (26) that only in-plane stresses ${\mathbf{\sigma }}_e$ can be obtained. The next order approximation is needed to obtain the out-of-plane stresses. The zeroth-order warping functions can be simply perturbed as

$${\mathit{w}}_{\left|\right|}={\overline{\mathit{v}}}_{\left|\right|},w_3={\overline{v}}_3+{\mathit{D}}_{\perp }\mathbf{\chi }$$

(29)

where $\mathbf{\chi }={\left[\mathbf{\epsilon }\mathbf{\kappa }\right]}^{\mathrm{T}},{\mathit{D}}_{\perp }=\left[-\frac{{\mathit{D}}_{\mathrm{et}}^{\mathrm{T}}}{{\mathit{D}}_t}-x_3\frac{{\mathit{D}}_{\mathrm{et}}^{\mathrm{T}}}{{\mathit{D}}_t}\right]$, ${\overline{\mathit{v}}}_{\left|\right|}$ and ${\overline{v}}_3$ denote in-plane and out-of-plane perturbed fluctuating functions, respectively.

The leading terms for the first-order approximation of variational statement can be obtained by substituting Equation (29) back into Equation (14) as

$$\begin{array}{c}2{\mathrm{\Pi }}_1=〈{\overline{\mathit{v}}}_{\left|\right|,3}^{\mathrm{T}}{\mathit{D}}_s{\overline{\mathit{v}}}_{\left|\right|,3}+{\mathit{D}}_t{\overline{v}}_{3,3}^2+2{\overline{\mathit{v}}}_{\left|\right|}^{\mathrm{T}}{\mathit{C}}_{\left|\right|,3}{\mathbf{\chi }}_{,\alpha }+2{\overline{\mathit{v}}}_{\left|\right|}^{\mathrm{T}}{\mathit{D}}_s{\partial }_t{\mathit{D}}_{\perp }{\mathbf{\chi }}_{,\alpha }-2{\overline{\mathit{v}}}_{\left|\right|}^{-\mathrm{T}}{\mathit{p}}_{\left|\right|}〉\hfill \\ -2{\mathit{v}}_{\left|\right|}^{-\mathrm{T}}{\mathbf{\tau }}_{\left|\right|}-2{\mathit{v}}_{\left|\right|}^{\mathrm{T}}{\mathbf{\beta }}_{\left|\right|}\hfill \end{array}$$

(30)

The Euler–Lagrange equation of Equation (30) is

$${\left({\mathit{D}}_s{\overline{\mathit{v}}}_{\left|\right|,3}+D_s{\partial }_t{\mathit{D}}_{\perp }{\mathbf{\chi }}_{,\alpha }\right)}_{,3}={\mathit{C}}_{\left|\right|,3}{\partial }_t{\mathit{D}}_{\perp }{\mathbf{\chi }}_{,\alpha }+{\mathit{g}}_{,3}+{\mathbf{\lambda }}_{\left|\right|}$$

(31)

where ${\mathit{C}}_{\left|\right|}=-{\partial }_e^{\mathrm{T}}\left[{\mathit{D}}_{\left|\right|}{x}_3{\mathit{D}}_{\left|\right|}\right],{\mathit{g}}_{,3}=-{\mathit{p}}_{\left|\right|}$.

It can be easily observed that ${\overline{v}}_3$ is decoupled from ${\overline{\mathit{v}}}_{\left|\right|}$ and only has a trivial solution. The solution of ${\overline{\mathit{v}}}_{\left|\right|}$ can be obtained from Equation (31) as

$${\overline{\mathit{v}}}_{\left|\right|}=\left({\overline{\mathit{C}}}_{\left|\right|}+{\mathit{L}}_{\alpha }\right){\mathbf{\chi }}_{,\alpha }+\overline{\mathit{g}}$$

(32)

where

$$\begin{array}{c}{\overline{\mathit{C}}}_{\left|\right|,3}={\mathit{D}}_s^{-1}{\mathit{C}}_{\left|\right|},〈{\overline{\mathit{C}}}_{\left|\right|}〉=0,{\overline{\mathit{g}}}_3={\mathit{D}}_s^{-1}\overline{\mathit{g}},\langle \overline{\mathit{g}}\rangle =0\hfill \\ {\mathit{L}}_{\alpha }{\mathbf{\chi }}_{,\alpha }=c_{\left|\right|}/h,{\overline{\mathit{C}}}_{\left|\right|}={\mathit{C}}_{\left|\right|}+\frac{x_3}{h}{\mathit{D}}_{\alpha }^{\mp }-\frac{1}{2}{\mathit{D}}_{\left|\right|}^{\pm }-{\mathit{D}}_s{\mathit{e}}_{\alpha }{\mathit{D}}_{\perp }\hfill \\ \overline{\mathit{g}}=\mathit{g}+\frac{x_3}{h}{\mathit{g}}^{\mp }-\frac{1}{2}{\mathit{g}}^{\pm }\hfill \end{array}$$

(33)

with ${\left(\right)}^{\pm }={\left(\right)}^{+}+{\left(\right)}^{-},{\left(\right)}^{\mp }={\left(\right)}^{-}-{\left(\right)}^{+}$.

So far, the asymptotic correction solution of the strain energy per unit area of the plate in the first-order approximation is

$$2{\mathrm{\Pi }}_1={\mathbf{\chi }}^{\mathrm{T}}{\overline{\mathit{D}}}_e\mathbf{\chi }+{\mathbf{\chi }}_{,\alpha }^{\mathrm{T}}{\mathit{B}}_{\alpha \beta }{\mathbf{\chi }}_{,\beta }-2{\mathbf{\chi }}^{\mathrm{T}}\mathit{F}$$

(34)

where

$$\begin{array}{c}{\mathit{B}}_{\alpha \beta }=〈D_{s_{\alpha \beta }}{\mathit{D}}_{\perp }^{\mathrm{T}}{\mathit{D}}_{\perp }-{\overline{\mathit{C}}}_{\alpha }^{\mathrm{T}}{\mathit{D}}_s^{-1}{\overline{\mathit{C}}}_{\beta }〉+{\mathit{L}}_{\alpha }{\left(\langle \overline{p}\rangle +〈{\mathit{p}}_{\left|\right|}〉\right)}_{,\alpha }\hfill \\ \mathit{F}=〈{\mathit{D}}_{\perp }^T{p}_3〉-〈{\overline{\mathit{C}}}_{\parallel }^T{\mathit{D}}_s^{-1}{\mathbf{g}}_s〉-{\mathit{L}}_{\alpha }{\left(\langle \overline{p}\rangle +〈{\mathit{p}}_{\parallel }〉\right)}_{,\alpha }\hfill \end{array}$$

(35)

The variational statement of Equation (34) governs the macroscopic behavior of the plate, and it only involves the 2-D field variables in terms of the macro-coordinates $x_1$ and $x_2$. Therefore, 2D-EPM can replace the original CCCP for the global analysis.

2.2.3. Recovery 3D Local Fields

The fidelity of the equivalent model should be evaluated based on how well it can predict the 3D local fields for the original 3D structure. Therefore, it is necessary to provide a recovery relationship to complete the equivalent model so that the results is comparable to those of the original 3D model.

The 3D displacement field can be recovered by using Equation (2) as

$$u_i=v_i+x_3\left(C_{3i}-{\delta }_{3i}\right)+C_{ji}{w}_j$$

(36)

where $u_i$ and $v_i$ are displacements of three-dimensional plate and two-dimensional equivalent plate, respectively, and $C_{ji}$ is the cosine component of the transformation matrix from macro coordinates to micro coordinates.

The 3D strain field can be recovered from Equation (8) as

$${\mathbf{\Gamma }}_e=\mathbf{\epsilon }+x_3\mathbf{\kappa },2{\mathbf{\Gamma }}_s={\overline{\mathit{v}}}_{\left|\right|3}+{\partial }_t{\mathit{D}}_{\perp ,\alpha },{\mathbf{\Gamma }}_t={\mathit{D}}_{\perp ,3}\mathbf{\chi }$$

(37)

At last, the 3D stress field can be recovered as

$$\begin{array}{c}{\mathbf{\sigma }}_e={\left[{\sigma }_{11}{\sigma }_{12}{\sigma }_{22}\right]}^{\mathrm{T}}={\mathit{D}}_{\left|\right|}\left(\mathbf{\epsilon }+x_3\mathbf{\kappa }\right)+{\mathit{D}}_e{\partial }_e{\overline{\mathit{v}}}_{\left|\right|,\alpha }\hfill \\ 2{\mathbf{\sigma }}_s={\left[{\sigma }_{13}{\sigma }_{23}\right]}^{\mathrm{T}}={\mathit{D}}_s\left({\overline{\mathit{v}}}_{\left|\right|,3}+{\partial }_t{\mathbf{\chi }}_{,\alpha }\right)\hfill \\ {\mathbf{\sigma }}_t={\sigma }_{33}={\mathit{D}}_{et}^{\mathrm{T}}{\partial }_e{\overline{\mathit{v}}}_{\left|\right|,\alpha }\hfill \end{array}$$

(38)

2.3. Free-Vibration Analysis of CCCP Using 2D-EPM

The elastic curved surface differential equation of 2D-EPM under lateral load $\tau$ is

$$D_{11}\frac{{\partial }^{4}v}{\partial x_1^4}+2D_{66}\frac{{\partial }^{4}v}{\partial x_1^2\partial x_2^2}+D_{22}\frac{{\partial }^{4}v}{\partial x_2^4}=\tau$$

(39)

It is assumed that the plate is in equilibrium under the lateral load $\tau$ at a certain moment, and the deflection of 2D-EPM at the equilibrium position is

$$v_a=v_a(x_1,x_2)$$

(40)

It is well known that 2D-EPM will be in a static state after moving for a period of time. Then, 2D-EPM will vibrate freely after the appropriate disturbing force is removed. The deflection of 2D-EPM at a certain moment under the condition of free vibration is

$$v_t=v_t(x_1,x_2,t)$$

(41)

In addition to the lateral load $\tau$, there is also the inertia force ${\tau }_j$ after the disturbing force is removed, and the equilibrium equation becomes

$$D_{11}\frac{{\partial }^4{v}_t}{\partial x_1^4}+2D_{66}\frac{{\partial }^4{v}_t}{\partial x_1^2\partial x_2^2}+D_{22}\frac{{\partial }^4{v}_t}{\partial x_2^4}=\tau +{\tau }_j$$

(42)

The acceleration of 2D-EPM is $\frac{{\partial }^{2}v}{\partial t^2}$. If the density and thickness of 2D-EPM are, respectively, ${\rho }_0$ and ${\delta }_0$, then the inertia force ${\tau }_j$ is

$${\tau }_j=-{\rho }_0{\delta }_0\frac{{\partial }^{2}v}{\partial t^2}$$

(43)

Substituting Equation (43) into Equation (42), and Equation (40) into Equation (39), the final differential equation can be obtained after subtracting Equation (39) from Equation (42) as

$$D_{11}\frac{{\partial }^4\left(v_t-v_a\right)}{\partial x_1^4}+2D_{66}\frac{{\partial }^4\left(v_t-v_a\right)}{\partial x_1^2\partial x_2^2}+D_{22}\frac{{\partial }^4\left(v_t-v_a\right)}{\partial x_2^4}=-{\rho }_0{\delta }_0\frac{{\partial }^2\left(v_t-v_a\right)}{\partial t^2}$$

(44)

The deflection of 2D-EPM satisfies the following conditions at any time:

$$v=v_t-v_a$$

(45)

Equation (44) can be simplified as

$$D_{11}\frac{{\partial }^{4}v}{\partial x_1^4}+2D_{66}\frac{{\partial }^{4}v}{\partial x_1^2\partial x_2^2}+D_{22}\frac{{\partial }^{4}v}{\partial x_2^4}+{\rho }_0{\delta }_0\frac{{\partial }^{2}v}{\partial t^2}=0$$

(46)

and its general solution is

$$v=\sum _{m=1}^{\infty }{v}_m=\sum _{m=1}^{\infty }\left(A_{m}cos\omega t+B_{m}sin\omega t\right)W(x,y)$$

(47)

where W is the free-vibration shape function and $\omega$ is the natural frequency (Hz).

As each natural frequency of 2D-EPM corresponds to its shape function, the differential equation can be obtained by substituting any of the natural frequency into Equation (46):

$$D_{11}\frac{{\partial }^{4}W}{\partial x_1^4}+2D_{66}\frac{{\partial }^{4}W}{\partial x_1^2\partial x_2^2}+D_{22}\frac{{\partial }^{4}W}{\partial x_2^4}=-{\omega }^2{\rho }_0{\delta }_{0}W$$

(48)

The nonzero solution satisfying W can be obtained as

$${\omega }^2=\frac{D_{11}\frac{{\partial }^{4}W}{\partial x_1^4}+D_{22}\frac{{\partial }^{4}W}{\partial x_2^4}+2D_{66}\frac{{\partial }^{4}W}{\partial x_1^2\partial x_2^2}}{{\rho }_0{\delta }_{0}W}$$

(49)

The free-vibration shape function of the four-side simple supporting plate can be defined as

$$W=sin\frac{m\pi x_1}{a}sin\frac{n\pi x_2}{b}$$

(50)

where m and n are half-wave numbers along $x_1$ and $x_2$ directions, respectively; a and b are the side lengths of 2D-EPM, respectively.

Substituting Equation (50) back into Equation (49), we can obtain

$${\omega }^2=\frac{D_{11}{\left(\frac{m\pi }{a}\right)}^4+D_{22}{\left(\frac{n\pi }{b}\right)}^4+2D_{66}{\left(\frac{m\pi }{a}\right)}^2{\left(\frac{n\pi }{b}\right)}^2}{{\rho }_0{\delta }_0}$$

(51)

Therefore, the natural frequencies of 2D-EPM can be obtained as

$$\omega ={\pi }^2\sqrt{\frac{D_{11}{\left(\frac{m}{a}\right)}^4+D_{22}{\left(\frac{n}{b}\right)}^4+2D_{66}{\left(\frac{m}{a}\right)}^2{\left(\frac{n}{b}\right)}^2}{{\rho }_0{\delta }_0}}$$

(52)

3. Model Validation

In this section, the accuracy and effectiveness of the equivalent model are verified by comparing with the results of static displacement, free-vibration, and global buckling analysis of the three-dimensional finite element model (3D-FEM).

The dimensions of unit cell are $l=30$ mm, $h=10$ mm, $t=1$ mm, $l_1=5$ mm, $l_2=5$ mm, and $l_3=10$ mm. The whole model of CCCP is the repetition of the unit cell in Figure 2b 20 times along $x_1$ and $x_2$ directions. The layup configuration of laminate is $[0/45/90/-45/0]\mathrm{s}$, and each ply thickness equals 0.1 mm. The material used is a unidirectional carbon fiber/epoxy resin composite (T300/7901), and the lamina properties are $E_{11}=E_{22}$ = 48.36 GPa, $E_{33}$ = 12.74 GPa, $G_{12}$ = 11.37 GPa, $G_{13}$ = $G_{23}$ = 4.59 GPa, $v_{12}$ = 0.078, $v_{13}$ = $v_{23}$ = 0.397, $\rho$ = 1670 $\mathrm{kg}/{\mathrm{m}}^3$. The effective plate properties of CCCP obtained by present model is shown in Figure 4.

3.1. Static Deformation Analysis

To investigate the static deformation of CCCP under different boundary conditions, the four cases shown in Figure 5 are considered.

Table 1. Comparison of static displacement between 3D-FEM and 2D-EPM under different boundary conditions (unit: mm).

Cases	Case 1	Case 2	Case 3	Case 4
3D-FEM
2D-EPM
Error %	1.03	4.28	2.05	5.21

lists the static displacements of 2D-EPM and 3D-FEM when a concentrated load of 100 N is applied to the center of the plate. It can be seen that the displacement distribution of 2D-EPM agrees with that of 3D-FEM under different boundary conditions. The boundary condition has a great influence on the relative errors. The smaller the boundary constraint is, the larger the displacement error is. However, the maximum error is less than 6% in Case 4, indicating the equivalent stiffness obtained by VAM is accurate, and the present 2D-EPM can well reflect the static displacement of CCCP under different boundary conditions.

3.2. Recovery 3D Field

As shown in Equations (36)–(38), the global displacements and strains of 2D-EPM should be imputed into the corresponding unit cell to recover the local field distribution. Note that most of the simplified models cannot predict the accurate local field distribution.

It can be seen from Figure 6 that the distribution of local stress within the unit cell is not uniform. The values of ${\sigma }_{11}$ and ${\sigma }_{22}$ are larger on the conical convexity, while other stresses are mainly concentrated at the intersection between the flat plate and the conical convexity. That is, these stresses on the conical convexity are very small, and the conical convexity cannot completely resist the load. It can be reasonably explained that the plate fails at the intersection of the flat plate and the conical convexity.

Figure 7 shows the distribution of the local fields within the unit cell closest to the midpoint of the plate in Case 1 (CCCC). It is seen that the displacement within the unit cell is symmetrically distributed, and the minimum displacement is located at the edge of the flat plate, while the maximum displacement is located at the intersection of the flat plate and the conical convexity. The changes of $U_1$ and $U_2$ are opposite on both sides of the conical convexity, and the maximum and minimum values of $U_1$ and $U_2$ are also located close to the intersection of the flat plate and the conical convexity.

Figure 8 shows the local stress and displacement distributions along Path 1 predicted by 2D-EPM and 3D-FEM. It can be seen that the local stress and displacement distributions predicted by the two models are in good agreement, and the maximum error is less than 2%. The location of stress concentration and the change trend can be obtained directly according to the local field distribution, which provides a reference for the damage analysis of the structure.

3.3. Global Buckling Analysis

In this section, the linear buckling behavior of CCCP is predicted by 2D-EPM and 3D-FEM under the assumption of small displacements. Four sides of CCCP are simple-supported and two opposite sides are pressed.

Table 2. Comparison of the first five buckling mode of CCCP in Case 3 predicted by different models (unit: N).

Orders	1	2	3	4	5
2D-EPM
3D-FEM
Error %	2.08	0.09	0.95	1.49	1.53

shows that the global buckling modes of 2D-EPM are consistent with those of 3D-FEM. For example, there are one, two, three, and four half-wave along the horizontal direction in the first, second, third, and fourth mode shape of 2D-EPM and 3D-FEM, respectively. The maximum error of critical buckling load is 2.08% in the first order buckling mode, indicating the equivalent model based on VAM has high accuracy in buckling analysis of CCCP.

3.4. Free-Vibration Analysis

In this section, the natural frequencies and vibration mode shapes of CCCP are predicted by 2D-EPM and 3D-FEM. The geometry and boundary conditions are the same as those in Section 3.3.

Table 3. Comparison of natural frequencies (Hz) and mode shapes of CCCP in Case 3 predicted by different models.

Orders	1	2	3	4	5
2D-EPM
3D-FEM
Error %	1.03	4.55	3.17	4.63	2.96

shows that the vibration modes predicted by 2D-EPM are consistent with those predicted by 3D-FEM, and the maximum error of natural frequency is within 5%. The accuracy of 2D-EPM based on VAM is verified from the aspect of vibration characteristics.

3.5. Comparison of Calculation Efficiency

To further demonstrate the advantages of the model, the computational efficiencies of the three-dimensional FE model and two-dimensional equivalent model are compared as shown in

Table 4. Comparison of computational efficiency between 3D-FEM and 2D-EPM.

Modelling Information		3D FEM	2D-EPM
			Unit Cell	2D Plate
Element type		C3D10	C3D4	S4R
Number of element		197,356	42,720	14,400
Number of nodes		394,157	59,827	14,641
Running time	Static analysis	5 min 30 s	10 s	30 s
	Buckling analysis	10 min 30 s		25 s
	Free-vibration analysis	7 min 15 s		25 s

. It can be observed that the present model is more time-efficient and cost-efficient than 3D-FEM in performing static, global buckling, and free-vibration analysis of CCCP. Note that the geometry of CCCP studied in this article is relatively simple, and the advantages of 2D-EPM can be reflected in more complex composite structures.

4. Parameter Analysis

4.1. Influence of Structural Parameters

The structural parameters mainly include $l_1,l_3,h$ and t as shown in Figure 3. In this section, the influence of structural parameters on the effective performance of CCCP is analyzed by using control variable method.

Figure 9a shows the influence of $l_1$ on the equivalent stiffness of CCCP. The parameters $l_3$ and h remain unchanged when considering the influence of $l_1$, so the total length of the unit cell is variable. It can be observed that the equivalent tensile stiffness of CCCP increases with increasing of $l_1$, while the equivalent bending stiffness gradually decreases nonlinearly. The main reason is that with the increase of the spacing between the conical convex part, the concave-convex characteristics of CCCP become less obvious, and the geometry of the structure is more similar to an ordinary flat plate.

Figure 9b shows that the first three buckling loads and natural frequencies of CCCP decrease with the increase of $l_1$. The critical buckling load of CCCP with $l_1$ = 2 mm is about three times that of $l_1$ = 10 mm, while the difference of natural frequency is about double, indicating that $l_1$ has a great influence on the buckling load and natural frequency of CCCP.

Figure 10a shows the influence of $l_3$ on the equivalent stiffness of CCCP. It can be observed that the width of convexity $l_3$ has little influence on the equivalent tensile and bending stiffness of CCCP. The curves of $A_{11}$ and $A_{22}$ decrease slightly, while the curve of $D_{11}$ shows an gentle and irregular change. The buckling loads and natural frequencies of CCCP shown in Figure 10b increase slightly with the increasing $l_3$. Therefore, the width of convexity $l_3$ has little influence on the effective performance of CCCP.

Figure 11a shows the influence of height h on the equivalent stiffness of CCCP. It can be observed that with the increase in convex height h, the equivalent tensile stiffness of CCCP decreases within a small range. The bending stiffness reaches the minimum value at h = 8 mm and then increases gradually with increasing of h. It can be seen form Figure 11b that the critical buckling load and natural frequency of CCCP increase to a certain extent with increasing of h, indicating the height of the convexity has a certain influence on the overall performance of the structure.

Figure 12a shows the influence of thickness t on the equivalent stiffness of CCCP. It can be seen that the equivalent tensile stiffness of CCCP increases linearly with the increasing thickness t, while the equivalent bending stiffness increases nonlinearly, which has a great influence on the buckling load and natural frequency of CCCP as shown in Figure 12b. The reason is that with the plate thickness t increases, the orthogonality of CCCP becomes non-obvious.

4.2. The Effect of Layup Configuration

4.2.1. The Influence of Ply Angles

In this section, the influence of ply angles on the effective performance of CCCP is analyzed. The layup configuration of laminate is set as ${[0/\theta /0/\theta /0]}_{\mathrm{s}}$, where the value of $\theta$ is from $0^{\circ }$ to ${90}^{\circ }$ with an interval of ${15}^{\circ }$.

Figure 13a shows the influence of ply angle on equivalent stiffness of CCCP. It can be observed that the curves of $A_{11}$ and $D_{11}$ show a nonlinear downward trend with the increase of the ply angles, while the curves of $A_{22}$ and $D_{22}$ show a small increase between ply angles of $0^{\circ }$ and ${45}^{\circ }$, and a nonlinear increase between ply angles of ${45}^{\circ }$ and ${90}^{\circ }$. It can be observed from Figure 13b that the influence of ply angle on buckling load is not obvious, while the second-order natural frequency reaches the maximum value when ply angles is ${90}^{\circ }$.

4.2.2. The Influence of the Proportion of $0^{\circ }$ Ply

In this section, CCCP with the combination of $0^{\circ }$ and ${45}^{\circ }$ ply is adopted to study the influence of the proportion of $0^{\circ }$ ply on the effective stiffness, buckling, and natural frequencies. It can be observed from Figure 14a that the values of $A_{11}$ and $D_{11}$ gradually increase with the increase in the proportion of $0^{\circ }$ ply, while other stiffness gradually decrease. The buckling load of CCCP in Figure 14b gradually increases with the increase in the proportion of $0^{\circ }$ ply, and reaches its maximum value when the proportion of $0^{\circ }$ ply is 80%. The first two natural frequencies of CCCP increase and then decrease with increasing of the proportion of $0^{\circ }$ ply. Therefore, the buckling resistance of CCCP can be improved by optimizing the layup configuration in engineering application, so that the effective plate properties of CCCP can be brought into full play.

5. Bidirectional CCCP

The structure introduced in this section is a bidirectional CCCP, and the convex direction along each row and column of the structure is alternately arranged up and down periodically as shown in Figure 15, which can be used to improve heat transfer effect and reduce pressure drop. Each convexity and its periphery is regarded as a small element, and only four adjacent elements (2 × 2) need to be studied due to the periodicity of bidirectional CCCP.

The geometrical sizes of unit cell are: h = 10 mm, t = 1.0 mm, $l_1$ = 5 mm, $l_2$ = 5 mm, $l_3$ = 10 mm, periodic length l = 60 mm. The material properties and boundary conditions are the same as those in Section 3.1. The equivalent stiffness matrix of the bidirectional CCCP obtained by the present model is shown in Figure 16 for reference.

Table 5. Comparison of displacements of bidirectional CCCP under different boundary conditions predicted by 2D-EPM and 3D-FEM.

Cases	Case 1	Case 2	Case 3	Case 4
3D-FEM
2D-EPM
Error %	1.41	3.10	0.69	0.18

shows the displacements of bidirectional CCCP predicted by 3D-FEM and 2D-EPM when the concentrated load of 100 N is applied at the center of the plate. It can be seen that the displacements of bidirectional CCCP predicted by 2D-EPM and 3D-FEM under four different boundary conditions are basically the same, and the maximum displacement deviation is only 3.10%. Therefore, the VAM-based equivalent model can be used to evaluate the static behavior of bidirectional CCCP with confidence.

Based on the obtained equivalent stiffness, the buckling critical eigenvalues of the equivalent model obtained by linear buckling analysis are compared with those of 3D-FEM.

Table 6. Comparison of the first five global buckling modes and critical loads (N) of bidirectional CCCP in Case 3 predicted by different models.

Orders	1	2	3	4	5
2D-EPM
3D-FEM
Error %	2.78	3.56	2.98	2.48	2.36

shows that the global buckling modes of 2D-EPM are consistent with those of 3D-FEM, and the buckling critical load error of each mode is less than 4%. Therefore, the VAM-based equivalent model has high accuracy in global buckling analysis of complex convex-concave plates.

To fully verify the correctness of the equivalent model of the bidirectional CCCP, the vibration modes and natural frequencies of 2D-EPM are analyzed by using the SFSF boundary conditions.

Table 7. Comparison of the first five vibration modes and natural frequencies (Hz) of bidirectional CCCP under the SFSF boundary conditions predicted by different models.

Orders	1	2	3	4	5
2D-EPM
3D-FEM
Error %	1.76	2.15	2.66	2.34	2.59

shows the comparison results of the first five vibration modes and natural frequencies between 3D-FEM and 2D-EPM of bidirectional CCCP. It can be seen that the vibration modes of 2D-EPM are consistent with those of 3D-FEM, and the natural frequency error of each mode is less than 2.6%, indicating the VAM-based equivalent model can be confidently used to predict the free vibration of complex convex-concave plates.

It can be concluded from the comparative analysis of Table 1, Table 2, Table 5 and Table 6 that compared with CCCP, the static displacement of bidirectional CCCP decreases and the buckling load increases under the same loading and boundary conditions, indicating the stiffness and stability of the bidirectional CCCP are greatly improved, which is consistent with the results of the previous equivalent stiffness analysis.

6. Conclusions

In this work, a VAM-based equivalent model (2D-EPM) is established to predict the effective performance of CCCP. As the approximate energy of 2D-EPM is as close as possible to that of the original 3D plate, it can be used to replace the original CCCP for the global analysis. The following conclusions can be obtained.

(1) The global displacements, buckling modes, and vibration modes of 2D-EPM under different boundary conditions are coincident with those of 3D-FEM, verifying the accuracy of the VAM-based 2D-EPM. Furthermore, the details of local field distribution within unit cell of CCCP are well captured. Compared with CCCP, the static displacement of bidirectional CCCP decreases and the buckling load are significantly increased, indicating the stiffness and stability of the bidirectional CCCP are greatly improved.

(2) The parameter analysis shows that the structural anisotropy of CCCP is more obvious with increasing of the convex height h, and the equivalent tensile stiffness decreases gradually, while the change of the equivalent bending stiffness is the opposite of the trend. The equivalent tensile stiffness increases gradually, and the equivalent bending stiffness decreases nonlinearly with increasing of the spacing between adjacent convexities, which may be because the anisotropy of CCCP is not obvious when the spacing between adjacent convexities increases.

(3) The ply angle of composite laminate has a great influence on the equivalent stiffness of CCCP. With increasing of the ply angle, the equivalent tensile and bending stiffness along the $x_1$ direction present a nonlinear decrease, while the equivalent tensile and bending stiffness along the $x_2$ direction gradually increases, but it has no significant influence on the buckling critical load and natural frequency of CCCP. With the increase in proportion of $0^{\circ }$ ply, the change of tensile and bending stiffness show a trend of linear increase, and the buckling critical load is also increased to a certain extent.

In short, compared with 3D-FEM, the present model has high accuracy for static, buckling, and free vibration analysis of CCCP. At the same time, the DOFs of the 2D-EPM is greatly reduced, resulting in a great improvement in computational efficiency. The present model is particularly useful in the early design stage of composite convex-concave structures where many design trade-offs need to be made over a vast design space in terms of material selection, ply angles, and geometric parameters.