Second-Order Multiscale Analysis and Algorithms for Composite Mindlin Plate and Timoshenko Beam

1. Introduction

The transition between precise material characterization and reliable structural application is a core challenge in the development of advanced composite materials. To achieve this goal, we need to establish a unified model that can accurately describe the complex mechanical behavior of composite materials under different structural forms. In previous works, the composite Kirchhoff plate and Euler beam theories have been widely applied in the thin structures [1,2]. However, due to the assumption that cross-section after deformation is perpendicular to the neutral axis and the lateral shear effect is ignored, the deviation is relatively large when simulating moderate thick composite materials. In this case, it is necessary to consider the composite Mindlin plate [3] and Timoshenko beam models, which can recover the Kirchhoff plate and Euler beam models as the thickness decreases. It has provided convenience for the efficient characterization and simulation of the effectiveness of modern composite materials.

Modern composite materials, composed of two or more different materials, such as ultra-high temperature ceramic composite plate reinforced with superalloys, carbon fiber reinforced epoxy resin composite plate [4,5], are widely used in aerospace, automotive and other high-performance engineering fields. These materials are valued for their outstanding rigidity and superior corrosion resistance, demonstrating remarkable adaptability to specialized application requirements. From a mathematical point of view, solving these equations requires dealing with rapidly varying coefficients. Due to the large scale of calculations, the conventional method incurs significant computational costs. Traditional formulations of plate and beam theory often fall short in capturing the complex behavior of structures with multiple length scales. A practical method proposed to simplify these calculations is the homogenization method [6]. By applying multi-scale expansion, homogenized solutions and cell functions can be obtained to describe both macroscopic and microscopic properties. This method improves the computational efficiency by homogenizing the microstructure features into effective macroscopic parameters, avoiding the expensive discretization at the microscale. Furthermore, the oscillating behavior of solutions can be effectively corrected. Based on the establishment of mathematical theories about the homogenization method [7,8,9], the multi-scale finite element method (MsFEM) [10], the heterogeneous multi-scale method (HMM) [11,12], and the variational multi-scale method (VMM) have been proposed and developed. Noticing that only first-order asymptotic analysis was established, Cui and Cao considered second-order asymptotic analysis [13]. Applying this method to simulate the deformation of composite materials is more accurate, and corresponding finite element algorithms have been established [14,15,16]. This SOTS asymptotic analysis is now applied to a variety of problems, such as thermoelasticity [17,18,19], elasticity problems [20,21], and so on.

However, significant technical difficulties arise in directly applying this asymptotic expansion to the classical theories of thin structures, namely the Kirchhoff plate and Euler beam problems. In the process of asymptotic expansion and analysis, these problems are described by fourth-order partial differential equations. After conducting an asymptotic analysis of models, it has been demonstrated that three-dimensional composite plate can be simplified to two-dimensional uniform Kirchhoff plate as the thickness h and $\epsilon$ approach zero [22]. Lewiński and Józef used the homogenization method and two-scale expansion method to solve the composite Mindlin plate problem in two dimensions [23]. Zhu and Cui provided a new Second-Order Two-Scale(SOTS) method for the Mindlin plate theory, and the laminated homogenization plate was analyzed [24]. Wang and Cui considered the bending behaviors of composite plate with 3-D periodic configuration and established a SOTS method, proving the approximation property of the solution by energy norm estimation [25]. Huang et al. established a two-scale asymptotic homogenization method for periodic composite Kirchhoff plate and Euler beam [26,27]. Numerical comparisons demonstrated that the current method is both physically acceptable and highly accurate.

As is known that directly establishing the SOTS asymptotic model for the fourth-order Kirchhoff plate and Euler beam directly is not easy, so it is reasonable to alternatively perform this high-order expansion on the second-order equations for Mindlin plate and Timoshenko beam first, and then proceed to consider approximating the corresponding composite Kirchhoff plate and Euler beam by letting the thickness h approach zero. We believe that this two-step strategy can not only obtain the asymptotic behavior of the Mindlin plate and Timoshenko beam, but also capture the local variations of the corresponding Kirchhoff plate and Euler beam. This is the principal idea and the main interest of our work, which are not seen elsewhere to the best of our knowledge. It should also be noted that when h is very small, the zero shear strain constraint leads to the shear locking [28,29,30], and we will discuss how to avoid this phenomenon in the asymptotic computing.

For these purposes mentioned above, the paper is structured as follows: The governing equations for the composite Mindlin plate and Timoshenko beam are presented in Section 2, and the corresponding SOTS asymptotic expansions are derived in Section 3. The finite element algorithm and symbol definitions for the second-order asymptotic plate and beam models are presented in Section 4. Several numerical examples of the SOTS asymptotic expansions for the Mindlin plate and the Timoshenko beam are given. The approaching behavior for the thick-thin theory of the plate and beam for the stationary problem are shown in Section 5. The conclusion remarks are given in Section 6. For simplicity, the Einstein summation convention is used, where the repeated subscript implies summation.

2. Governing Equations

The configuration of the composite plate with thickness h is typically shown in Figure 1(a). Different constituent materials are periodically distributed in the $x_1$ and $x_2$ directions and are homogeneous in the thickness direction. The normalized periodic cell is depicted in Figure 1(b).

Define the two-dimensional middle plane as $\mathrm{\Omega }$ and the unit cell as Y. With the micro and macro coordinates being $\mathbf{x}=[x_1,x_2]\in \mathrm{\Omega }$ and $\mathbf{y}=[y_1,y_2]\in Y$, respectively, a small parameter $\epsilon$ is defined to relate these variables as

$$\mathbf{y}={\displaystyle \frac{\mathbf{x}}{\epsilon }}-⌊{\displaystyle \frac{\mathbf{x}}{\epsilon }}⌋,$$

(1)

where $⌊·⌋$ is floor operator.

The governing equations for the composite Mindlin plate can be formulated as

$$\left\{\begin{array}{c}{\displaystyle -\frac{\partial V_1^{\epsilon }}{\partial x_1}-\frac{\partial V_2^{\epsilon }}{\partial x_2}=q_{\mathrm{M}},}\hfill \\ {\displaystyle \frac{\partial M_1^{\epsilon }}{\partial x_1}+\frac{\partial M_{12}^{\epsilon }}{\partial x_2}-V_1^{\epsilon }=0,}\hfill \\ {\displaystyle \frac{\partial M_{12}^{\epsilon }}{\partial x_1}+\frac{\partial M_2^{\epsilon }}{\partial x_2}-V_2^{\epsilon }=0,}\hfill \end{array}\right.$$

(2)

where $V_1^{\epsilon }$ and $V_2^{\epsilon }$ are the shear stresses, $M_1^{\epsilon }$ and $M_2^{\epsilon }$ are the bending moments, $M_{12}^{\epsilon }$ is the torsional moment, and $q_{\mathrm{M}}$ represents the external load on the plate.

The shear stress $V_1^{\epsilon }$ and $V_2^{\epsilon }$ can be expressed by

$$\left[\begin{array}{c}{V}_1^{\epsilon }\\ V_2^{\epsilon }\end{array}\right]={\displaystyle \frac{E^{\epsilon }{k}_{s}h}{2(1+{\nu }^{\epsilon })}}\left[\begin{array}{c}{\displaystyle \frac{\partial w^{\epsilon }}{\partial x_1}-{\beta }_1^{\epsilon }}\\ {\displaystyle \frac{\partial w^{\epsilon }}{\partial x_2}-{\beta }_2^{\epsilon }}\end{array}\right],$$

(3)

where $w^{\epsilon }$ is the deflection, ${\beta }_1^{\epsilon }$ and ${\beta }_2^{\epsilon }$ are the rotations of the cross section, $E^{\epsilon }$ is the Young’s modulus, $k_s$ is the shear correction factor, and ${\nu }^{\epsilon }$ is the Poisson’s ratio. $E^{\epsilon }$ and ${\nu }^{\epsilon }$ are periodic due to the configuration and can be expressed as

$$\begin{array}{c}\hfill {\displaystyle E^{\epsilon }\left(\mathbf{x}\right)=E\left(\frac{\mathbf{x}}{\epsilon }\right)=E\left(\mathbf{y}\right),{\nu }^{\epsilon }\left(\mathbf{x}\right)=\nu \left(\frac{\mathbf{x}}{\epsilon }\right)=\nu \left(\mathbf{y}\right),}\end{array}$$

(4)

where $E\left(\mathbf{y}\right)$ and $\nu \left(\mathbf{y}\right)$ are periodic functions defined in Y. The constitutive relation between bending moment and curvature is given by

$$\left[\begin{array}{c}{M}_1^{\epsilon }\\ M_2^{\epsilon }\\ M_{12}^{\epsilon }\end{array}\right]={\displaystyle \frac{E^{\epsilon }{h}^3}{12\left(1-{\left({\nu }^{\epsilon }\right)}^2\right)}}\left[\begin{array}{ccc}1& {\nu }^{\epsilon }& 0\\ {\nu }^{\epsilon }& 1& 0\\ 0& 0& {\displaystyle \frac{1-{\nu }^{\epsilon }}{2}}\end{array}\right]\left[\begin{array}{c}{\displaystyle -\frac{\partial {\beta }_1^{\epsilon }}{\partial x_1}}\\ {\displaystyle -\frac{\partial {\beta }_2^{\epsilon }}{\partial x_2}}\\ {\displaystyle -\frac{\partial {\beta }_1^{\epsilon }}{\partial x_2}-\frac{\partial {\beta }_2^{\epsilon }}{\partial x_1}}\end{array}\right].$$

(5)

For simplicity, clamped boundary conditions are imposed on $\partial \mathrm{\Omega }$ and the Mindlin plate problem is given as

$$\left\{\begin{array}{c}\hfill \begin{array}{cc}\hfill {\displaystyle -\frac{\partial }{\partial x_i}\left(k_s{G}^{\epsilon }h\left(\frac{\partial w^{\epsilon }}{\partial x_i}-{\beta }_i^{\epsilon }\right)\right)}& =q_{\mathrm{M}}\mathrm{in}\mathrm{\Omega },\hfill \\ \hfill {\displaystyle -\frac{\partial }{\partial x_j}\left(C_{ijkl}^{\epsilon }\frac{\partial {\beta }_k^{\epsilon }}{\partial x_l}\right)-k_s{G}^{\epsilon }h\left(\frac{\partial w^{\epsilon }}{\partial x_i}-{\beta }_i^{\epsilon }\right)}& =0\mathrm{in}\mathrm{\Omega },\hfill \\ \hfill w^{\epsilon }& =0\mathrm{on}\partial \mathrm{\Omega },\hfill \\ \hfill {\beta }_i^{\epsilon }& =0\mathrm{on}\partial \mathrm{\Omega },\hfill \end{array}\end{array}\right.$$

(6)

where $C_{ijkl}^{\epsilon }$ is defined by

$$C_{ijkl}^{\epsilon }={\displaystyle \frac{E^{\epsilon }{h}^3{\nu }^{\epsilon }}{12(1-{\left({\nu }^{\epsilon }\right)}^2)}}{\delta }_{ij}{\delta }_{kl}+{\displaystyle \frac{E^{\epsilon }{h}^3}{24\left(1+{\nu }^{\epsilon }\right)}}\left({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk}\right),$$

(7)

and $G^{\epsilon }$ is the shear modulus, related to $E^{\epsilon }$ and ${\nu }^{\epsilon }$ through $G^{\epsilon }={\displaystyle \frac{E^{\epsilon }}{2(1+{\nu }^{\epsilon })}}$.

The Timoshenko beam problem is equivalent to the dimensional reduction of the plate problem (6). Therefore, the composite Timoshenko beam problem takes the following form:

$$\left\{\begin{array}{c}\hfill \begin{array}{cc}\hfill {\displaystyle -\frac{d}{d{x}_1}\left(k_s{G}^{\epsilon }A\left(\frac{d{w}^{\epsilon }}{d{x}_1}-{\beta }^{\epsilon }\right)\right)}& =q_{\mathrm{T}}\mathrm{in}\mathrm{\Omega },\hfill \\ \hfill {\displaystyle -\frac{d}{d{x}_1}\left(E^{\epsilon }I\frac{d{\beta }^{\epsilon }}{d{x}_1}\right)-k_s{G}^{\epsilon }A\left(\frac{d{w}^{\epsilon }}{d{x}_1}-{\beta }^{\epsilon }\right)}& =0\mathrm{in}\mathrm{\Omega },\hfill \\ \hfill w^{\epsilon }& =0\mathrm{on}\partial \mathrm{\Omega },\hfill \\ \hfill {\beta }^{\epsilon }& =0\mathrm{on}\partial \mathrm{\Omega },\hfill \end{array}\end{array}\right.$$

(8)

where $w^{\epsilon }$ is the one-dimensional deflection, ${\beta }^{\epsilon }$ is the rotation, A is the sectional area, given by $A=b\times h$ with b being the width, and $q_{\mathrm{T}}$ is the external load on the beam.

3. The SOTS Asymptotic Analysis

First, we perform the SOTS asymptotic analysis for the composite Mindlin plate problems (6). Following the standard procedure, the deflection $w^{\epsilon }$, rotation ${\beta }_i^{\epsilon }$ can be expanded as

$$\begin{array}{cc}\hfill w^{\epsilon }\left(\mathbf{x}\right)& =w_0\left(\mathbf{x}\right)+\epsilon w_1(\mathbf{x},\mathbf{y})+{\epsilon }^2{w}_2(\mathbf{x},\mathbf{y})+O\left({\epsilon }^3\right),\hfill \\ \hfill {\beta }_i^{\epsilon }\left(\mathbf{x}\right)& ={\beta }_i^0\left(\mathbf{x}\right)+\epsilon {\beta }_i^1(\mathbf{x},\mathbf{y})+{\epsilon }^2{\beta }_i^2(\mathbf{x},\mathbf{y})+O\left({\epsilon }^3\right),i=1,2,\hfill \end{array}$$

(9)

with the homogenized solutions $w_0\left(\mathbf{x}\right)$, ${\beta }_i^0\left(\mathbf{x}\right)$, appended by the first-order correctors $w_1(\mathbf{x},\mathbf{y})$, ${\beta }_i^1(\mathbf{x},\mathbf{y})$ and the second-order terms $w_2(\mathbf{x},\mathbf{y})$, ${\beta }_i^2(\mathbf{x},\mathbf{y})$. Transforming the partial derivatives as

$${\displaystyle \frac{\partial }{\partial x_i}}\to {\displaystyle \frac{\partial }{\partial x_i}}+{\displaystyle \frac{1}{\epsilon }}{\displaystyle \frac{\partial }{\partial y_i}},$$

(10)

substituting Eqs. (9) into governing equations (6), and equating the power coefficients of $\epsilon$, the equations for the first three expansion terms are obtained as follows:

$$O\left({\epsilon }^{-1}\right):\left\{\begin{array}{c}{\displaystyle -\frac{\partial }{\partial y_i}\left(k_{s}Gh\frac{\partial w_1}{\partial y_i}\right)-\frac{\partial }{\partial y_i}\left(k_{s}Gh\left(\frac{\partial w_0}{\partial x_i}-{\beta }_i^0\right)\right)=0,}\hfill \\ {\displaystyle -\frac{\partial }{\partial y_j}\left(C_{ijkl}\frac{\partial {\beta }_k^1}{\partial y_l}\right)-\frac{\partial }{\partial y_j}\left(C_{ijkl}\frac{\partial {\beta }_k^0}{\partial x_l}\right)=0,}\hfill \end{array}\right.$$

(11)

$$O\left({\epsilon }^0\right):\left\{\begin{array}{c}{\displaystyle -\frac{\partial }{\partial x_i}\left(k_{s}Gh\left(\frac{\partial w_1}{\partial y_i}+\frac{\partial w_0}{\partial x_i}-{\beta }_i^0\right)\right)+\frac{\partial }{\partial y_i}\left(k_{s}Gh{\beta }_i^1\right)}\hfill \\ {\displaystyle -\frac{\partial }{\partial y_i}\left(k_{s}Gh\frac{\partial w_2}{\partial y_i}\right)-\frac{\partial }{\partial y_i}\left(k_{s}Gh\frac{\partial w_1}{\partial x_i}\right)=q_{\mathrm{M}},}\hfill \\ {\displaystyle -\frac{\partial }{\partial y_j}\left(C_{ijkl}\frac{\partial {\beta }_k^2}{\partial y_l}\right)-\frac{\partial }{\partial x_j}\left(C_{ijkl}\frac{\partial {\beta }_k^1}{\partial y_l}\right)-\frac{\partial }{\partial y_j}\left(C_{ijkl}\frac{\partial {\beta }_k^1}{\partial x_l}\right)}\hfill \\ {\displaystyle -\frac{\partial }{\partial x_j}\left(C_{ijkl}\frac{\partial {\beta }_k^0}{\partial x_l}\right)-k_{s}Gh\left(\frac{\partial w_0}{\partial x_i}+\frac{\partial w_1}{\partial y_i}-{\beta }_i^0\right)=0.}\hfill \end{array}\right.$$

(12)

From (11), the first-order correctors $w_1$ and ${\beta }_k^1$ are defined as

$$\begin{array}{c}\hfill {\displaystyle w_1(\mathbf{x},\mathbf{y})=N_{{\alpha }_1}\left(\mathbf{y}\right)\left({\displaystyle \frac{\partial w_0\left(\mathbf{x}\right)}{\partial x_{{\alpha }_1}}-{\beta }_{{\alpha }_1}^0\left(\mathbf{x}\right)}\right),{\beta }_k^1(\mathbf{x},\mathbf{y})=N_{{\alpha }_{1}km}\left(\mathbf{y}\right)\frac{\partial {\beta }_m^0\left(\mathbf{x}\right)}{\partial x_{{\alpha }_1}},}\end{array}$$

(13)

where the first-order cell functions $N_{{\alpha }_1}$ and $N_{{\alpha }_{1}km}$ are solutions to the following boundary value problems:

$$\left\{\begin{array}{c}{\displaystyle -\frac{\partial }{\partial y_i}\left(k_{s}Gh\frac{\partial N_{{\alpha }_1}}{\partial y_i}\right)=\frac{\partial }{\partial y_{{\alpha }_1}}\left(k_{s}Gh\right)\mathrm{in}Y,}\hfill \\ {\displaystyle {\int }_Y{N}_{{\alpha }_1}\left(\mathbf{y}\right)d\mathbf{y}=0,}\hfill \end{array}\right.$$

(14)

$$\left\{\begin{array}{c}{\displaystyle -\frac{\partial }{\partial y_j}\left(C_{ijkl}\frac{\partial N_{{\alpha }_{1}km}}{\partial y_l}\right)=\frac{\partial }{\partial y_j}\left(C_{ijm{\alpha }_1}\right)\mathrm{in}Y,}\hfill \\ {\displaystyle {\int }_Y{N}_{{\alpha }_{1}km}\left(\mathbf{y}\right)d\mathbf{y}=0.}\hfill \end{array}\right.$$

(15)

Next, taking the integral average over the unit cell Y in (12), and using the newly defined expressions for $w_1$ and ${\beta }_k^1$ from (13), the homogenized Mindlin plate problem can be defined as follows:

$$\left\{\begin{array}{c}\hfill \begin{array}{cc}\hfill {\displaystyle -\frac{\partial }{\partial x_i}\left(k_s{G}_{ij}^{0}h\left(\frac{\partial w_0}{\partial x_j}-{\beta }_j^0\right)\right)}& =q_{\mathrm{M}}^0\mathrm{in}\mathrm{\Omega },\hfill \\ \hfill {\displaystyle -\frac{\partial }{\partial x_j}\left(C_{ijkl}^0\frac{\partial {\beta }_k^0}{\partial x_l}\right)-k_s{G}_{ij}^{0}h\left(\frac{\partial w_0}{\partial x_j}-{\beta }_j^0\right)}& =0\mathrm{in}\mathrm{\Omega },\hfill \\ \hfill w_0& =0\mathrm{on}\partial \mathrm{\Omega },\hfill \\ \hfill {\beta }_i^0& =0\mathrm{on}\partial \mathrm{\Omega },\hfill \end{array}\end{array}\right.$$

(16)

where the homogenized coefficients $G_{ij}^0$ and $C_{ijkl}^0$, and the homogenized load $q_{\mathrm{M}}^0$ are given by

$$\begin{array}{cc}\hfill G_{ij}^0& ={\displaystyle \frac{1}{\left|Y\right|}}{\int }_{Y}G\left({\delta }_{ij}+{\displaystyle \frac{\partial N_j}{\partial y_i}}\right)d\mathbf{y},C_{ijkl}^0={\displaystyle \frac{1}{\left|Y\right|}}{\int }_Y\left(C_{ijkl}+C_{ijmn}{\displaystyle \frac{\partial N_{lmk}}{\partial y_n}}\right)d\mathbf{y},\hfill \\ \hfill q_{\mathrm{M}}^0& ={\displaystyle \frac{1}{\left|Y\right|}}{\int }_Y{q}_{\mathrm{M}}\mathrm{d}\mathbf{y},\hfill \end{array}$$

(17)

respectively, where $\left|Y\right|$ is the measure of the domain Y.

A high-order expansion is necessary to capture the local oscillation behavior of deflection and rotation. By substituting the homogenized equation (16) and the first-order expansion (13) into (12), the second-order correctors $w_2$ and ${\beta }_k^2$ have the following expressions:

$$\begin{array}{cc}\hfill w_2(\mathbf{x},\mathbf{y})& =N_{{\alpha }_1{\alpha }_2}\left(\mathbf{y}\right)\left({\displaystyle \frac{{\partial }^2{w}_0\left(\mathbf{x}\right)}{\partial x_{{\alpha }_1}\partial x_{{\alpha }_2}}}-{\displaystyle \frac{\partial {\beta }_{{\alpha }_1}^0\left(\mathbf{x}\right)}{\partial x_{{\alpha }_2}}}\right)-M_{{\alpha }_1{\alpha }_2}\left(\mathbf{y}\right){\displaystyle \frac{\partial {\beta }_{{\alpha }_2}^0\left(\mathbf{x}\right)}{\partial x_{{\alpha }_1}}},\hfill \\ \hfill {\beta }_k^2(\mathbf{x},\mathbf{y})& =N_{{\alpha }_1{\alpha }_{2}km}\left(\mathbf{y}\right){\displaystyle \frac{{\partial }^2{\beta }_m^0\left(\mathbf{x}\right)}{\partial x_{{\alpha }_1}\partial x_{{\alpha }_2}}}+L_{{\alpha }_{1}k}\left(\mathbf{y}\right)\left({\displaystyle \frac{\partial w_0\left(\mathbf{x}\right)}{\partial x_{{\alpha }_1}}}-{\beta }_{{\alpha }_1}^0\left(\mathbf{x}\right)\right),\hfill \end{array}$$

(18)

where the second-order cell functions $N_{{\alpha }_1{\alpha }_2}$, $M_{{\alpha }_1{\alpha }_2}$, $N_{{\alpha }_1{\alpha }_{2}km}$, and $L_{{\alpha }_{1}k}$ are solutions to the following cell problems:

$$\left\{\begin{array}{cc}\hfill -{\displaystyle \frac{\partial }{\partial y_i}}\left(k_{s}Gh{\displaystyle \frac{\partial N_{{\alpha }_1{\alpha }_2}}{\partial y_i}}\right)=& {\displaystyle \frac{\partial }{\partial y_{{\alpha }_2}}}\left(k_{s}Gh{N}_{{\alpha }_1}\right)+k_{s}Gh{\displaystyle \frac{\partial N_{{\alpha }_1}}{\partial y_{{\alpha }_2}}}+k_{s}Gh{\delta }_{{\alpha }_1{\alpha }_2}-k_s{G}_{{\alpha }_1{\alpha }_2}^{0}h{\displaystyle \frac{q_{\mathrm{M}}}{q_{\mathrm{M}}^0}}\mathrm{in}Y,\hfill \\ \hfill {\int }_Y{N}_{{\alpha }_1{\alpha }_2}\left(\mathbf{y}\right)d\mathbf{y}& =0,\hfill \end{array}\right.$$

$$\left\{\begin{array}{cc}\hfill -{\displaystyle \frac{\partial }{\partial y_i}}\left(k_{s}Gh{\displaystyle \frac{\partial M_{{\alpha }_1{\alpha }_2}}{\partial y_i}}\right)& ={\displaystyle \frac{\partial }{\partial y_i}}\left(k_{s}Gh{N}_{{\alpha }_{1}i{\alpha }_2}\right)\mathrm{in}Y,\hfill \\ \hfill {\int }_Y{M}_{{\alpha }_1{\alpha }_2}\left(\mathbf{y}\right)d\mathbf{y}& =0,\hfill \end{array}\right.$$

$$\left\{\begin{array}{cc}\hfill -{\displaystyle \frac{\partial }{\partial y_j}}\left(C_{ijkl}{\displaystyle \frac{\partial N_{{\alpha }_1{\alpha }_{2}km}}{\partial y_l}}\right)& ={\displaystyle \frac{\partial }{\partial y_j}}\left(C_{ijk{\alpha }_2}{N}_{{\alpha }_{1}km}\right)+C_{i{\alpha }_{2}kl}{\displaystyle \frac{\partial N_{{\alpha }_{1}km}}{\partial y_l}}+C_{i{\alpha }_{2}m{\alpha }_1}-C_{i{\alpha }_{2}m{\alpha }_1}^0\mathrm{in}Y,\hfill \\ \hfill {\int }_Y{N}_{{\alpha }_1{\alpha }_{2}km}\left(\mathbf{y}\right)d\mathbf{y}& =0,\hfill \end{array}\right.$$

$$\left\{\begin{array}{cc}\hfill -{\displaystyle \frac{\partial }{\partial y_j}}\left(C_{ijkl}{\displaystyle \frac{\partial L_{{\alpha }_{1}k}}{\partial y_l}}\right)& =k_{s}Gh{\delta }_{i{\alpha }_1}+k_{s}Gh{\displaystyle \frac{\partial N_{{\alpha }_1}}{\partial y_i}}-k_s{G}_{i{\alpha }_1}^{0}h\mathrm{in}Y,\hfill \\ \hfill {\int }_Y{L}_{{\alpha }_{1}k}\left(\mathbf{y}\right)d\mathbf{y}& =0,\hfill \end{array}\right.$$

Next, we move to the Timoshenko beam problem (8). Based on the results from the Mindlin plate problem, the expansions of $w^{\epsilon }$ and ${\beta }^{\epsilon }$ can be directly obtained as

$$\begin{array}{cc}\hfill w^{\epsilon }\left(x_1\right)=& w_0\left(x_1\right)+\epsilon w_1(x_1,y_1)+{\epsilon }^2{w}_2(x_1,y_1)+O\left({\epsilon }^3\right),\hfill \\ \hfill =& w_0\left(x_1\right)+\epsilon N_1\left(y_1\right)\left({\displaystyle \frac{d{w}_0\left(x_1\right)}{d{x}_1}}-{\beta }^0\left(x_1\right)\right)+{\epsilon }^2\left(N_{11}\left(y_1\right){\displaystyle \frac{d^2{w}_0\left(x_1\right)}{d{x}_1^2}}+M\left(y_1\right){\displaystyle \frac{d{\beta }^0\left(x_1\right)}{d{x}_1}}\right)\hfill \\ \hfill & +O\left({\epsilon }^3\right),\hfill \\ \hfill {\beta }^{\epsilon }\left(x_1\right)=& {\beta }^0\left(x_1\right)+\epsilon {\beta }^1(x_1,y_1)+{\epsilon }^2{\beta }^2(x_1,y_1)+O\left({\epsilon }^3\right),\hfill \\ \hfill =& {\beta }^0\left(x_1\right)+\epsilon H_1\left(y_1\right){\displaystyle \frac{d{\beta }^0\left(x_1\right)}{d{x}_1}}+{\epsilon }^2\left(H_{11}\left(y_1\right)\left({\beta }^0\left(x_1\right)-{\displaystyle \frac{d{w}_0\left(x_1\right)}{d{x}_1}}\right)+M_{11}\left(y_1\right){\displaystyle \frac{d^2{\beta }^0\left(x_1\right)}{d{x}_1^2}}\right)\hfill \\ \hfill & +O\left({\epsilon }^3\right),\hfill \end{array}$$

(19)

where the first-order cell functions $N_1\left(y_1\right)$ and $H_1\left(y_1\right)$ respectively satisfy the following problems:

$$\left\{\begin{array}{c}{\displaystyle -\frac{d}{d{y}_1}\left(k_{s}GA\frac{d{N}_1}{d{y}_1}\right)=\frac{d}{d{y}_1}\left(k_{s}GA\right)\mathrm{in}Y,}\hfill \\ {\displaystyle {\int }_Y{N}_1\left(y_1\right)d{y}_1=0,}\hfill \end{array}\right.$$

(20)

$$\left\{\begin{array}{c}{\displaystyle -\frac{d}{d{y}_1}\left(EI\frac{d{H}_1}{d{y}_1}\right)=\frac{d}{d{y}_1}\left(EI\right)\mathrm{in}Y,}\hfill \\ {\displaystyle {\int }_Y{H}_1\left(y_1\right)d{y}_1=0,}\hfill \end{array}\right.$$

(21)

The homogenized solutions $w_0$ and ${\beta }^0$ are obtained from

$$\left\{\begin{array}{c}\hfill \begin{array}{cc}\hfill {\displaystyle -\frac{d}{d{x}_1}\left(k_s{G}^{0}A\left(\frac{d{w}_0}{d{x}_1}-{\beta }^0\right)\right)}& =q_{\mathrm{T}}^0\mathrm{in}\mathrm{\Omega },\hfill \\ \hfill {\displaystyle -\frac{d}{d{x}_1}\left(E^{0}I\frac{d{\beta }^0}{d{x}_1}\right)-k_s{G}^{0}A\left(\frac{d{w}_0}{d{x}_1}-{\beta }^0\right)}& =0\mathrm{in}\mathrm{\Omega },\hfill \\ \hfill w_0& =0\mathrm{on}\partial \mathrm{\Omega },\hfill \\ \hfill {\beta }^0& =0\mathrm{on}\partial \mathrm{\Omega },\hfill \end{array}\end{array}\right.$$

(22)

where the homogenized coefficients $G^0$ and $E^0$, and $q_{\mathrm{T}}^0$ from (22) are respectively defined by

$$G^0={\displaystyle \frac{1}{\left|Y\right|}}{\int }_{Y}G\left(1+{\displaystyle \frac{d{N}_1}{d{y}_1}}\right)d{y}_1,E^0={\displaystyle \frac{1}{\left|Y\right|}}{\int }_{Y}E\left(1+{\displaystyle \frac{d{H}_1}{d{y}_1}}\right)d{y}_1,q_{\mathrm{T}}^0={\displaystyle \frac{1}{\left|Y\right|}}{\int }_Y{q}_{\mathrm{T}}d{y}_1.$$

(23)

The second-order cell functions $N_{11}\left(y_1\right)$, $M\left(y_1\right)$, $H_{11}\left(y_1\right)$, and $M_{11}\left(y_1\right)$ are solutions to the problems:

$$\begin{array}{cc}\hfill & \left\{\begin{array}{cc}\hfill & {\displaystyle -\frac{d}{d{y}_1}\left(k_{s}GA\frac{d{N}_{11}}{d{y}_1}\right)=k_{s}GA\frac{d{N}_1}{d{y}_1}+\frac{d}{d{y}_1}\left(k_{s}GA{N}_1\right)+k_{s}GA-k_s{G}^{0}A\frac{q_{\mathrm{T}}}{q_{\mathrm{T}}^0}\mathrm{in}Y,}\hfill \\ \hfill & {\int }_Y{N}_{11}\left(y_1\right)d{y}_1=0,\hfill \end{array}\right.\hfill \\ \hfill & \left\{\begin{array}{cc}\hfill & {\displaystyle -\frac{d}{d{y}_1}\left(k_{s}GA\frac{dM}{d{y}_1}\right)=\frac{d}{d{y}_1}\left(k_{s}GA{N}_1\right)+\frac{d}{d{y}_1}\left(k_{s}GA{H}_1\right)\mathrm{in}Y,}\hfill \\ \hfill & {\int }_{Y}M\left(y_1\right)d{y}_1=0,\hfill \end{array}\right.\hfill \\ \hfill & \left\{\begin{array}{cc}\hfill & -{\displaystyle \frac{d}{d{y}_1}}\left(EI{\displaystyle \frac{d{H}_{11}}{d{y}_1}}\right)=-k_{s}GA\left({\displaystyle \frac{d{N}_1}{d{y}_1}}+1\right)+k_s{G}^{0}A\mathrm{in}Y,\hfill \\ \hfill & {\int }_Y{H}_{11}\left(y_1\right)d{y}_1=0,\hfill \end{array}\right.\hfill \\ \hfill & \left\{\begin{array}{cc}\hfill & -{\displaystyle \frac{d}{d{y}_1}}\left(EI{\displaystyle \frac{d{M}_{11}}{d{y}_1}}\right)=EI{\displaystyle \frac{d{H}_1}{d{y}_1}}+{\displaystyle \frac{d}{d{y}_1}}\left(EI{H}_1\right)+EI-E^{0}I\mathrm{in}Y,\hfill \\ \hfill & {\int }_Y{M}_{11}\left(y_1\right)d{y}_1=0.\hfill \end{array}\right.\hfill \end{array}$$

4. Multi-Scale Finite Element Algorithm

According to the expansions of the deflections and rotations, we can establish the finite element(FE) algorithms for the plate and beam problems. The multiscale computational flowchart for Mindlin plate is shown in Figure 2. Two kinds of meshes are generated before the multiscale computations: cell mesh and homogenized mesh. All the cell functions are solved on the cell mesh, and the homogenized solutions are obtained on the homogenized mesh. First, let us define the first-order and second-order approximate solutions as follows:

$$\begin{array}{cc}\hfill & w^{\epsilon ,1}=w_0+\epsilon w_1,w^{\epsilon ,2}=w^{\epsilon ,1}+{\epsilon }^2{w}_2,\hfill \\ \hfill & {\beta }_i^{\epsilon ,1}={\beta }_i^0+\epsilon {\beta }_i^1,{\beta }_i^{\epsilon ,2}={\beta }_i^{\epsilon ,1}+{\epsilon }^2{\beta }_i^2,\hfill \end{array}$$

(24)

with the homogenized shear stresses $V_i^0$ and the first-order and second-order approximate shear stresses $V_i^{\epsilon ,1}$ and $V_i^{\epsilon ,2}$, respectively, given by

$$\begin{array}{cc}\hfill & V_i^0={\displaystyle \frac{E^{\epsilon }{k}_{s}h}{2(1+{\nu }^{\epsilon })}}\left({\displaystyle \frac{\partial w_0}{\partial x_i}}-{\beta }_0\right),\hfill \\ \hfill & V_i^{\epsilon ,1}=V_i^0+\epsilon {\displaystyle \frac{E^{\epsilon }{k}_{s}h}{2(1+{\nu }^{\epsilon })}}\left({\displaystyle \frac{\partial w_1}{\partial x_i}}-{\beta }_i^1\right),\hfill \\ \hfill & V_i^{\epsilon ,2}=V_i^{\epsilon ,1}+{\epsilon }^2{\displaystyle \frac{E^{\epsilon }{k}_{s}h}{2(1+{\nu }^{\epsilon })}}\left({\displaystyle \frac{\partial w_2}{\partial x_i}}-{\beta }_i^2\right).\hfill \end{array}$$

(25)

The FE procedures for plate problem are illustrated as follows:

1. Preprocessing: Set the material coefficients $E^{\epsilon }$ and ${\nu }^{\epsilon }$, define the periodic cell Y and the homogenized domain $\mathrm{\Omega }$, perform the spatial discretization and generate the computational meshes. To illustrate the efficiency of the SOTS model proposed in this paper, we also partition a sufficiently refined mesh and conduct the classical computations for $w^{\epsilon }$ and ${\beta }_i^{\epsilon }$ as the reference solution for comparison.

2. Computations of the cell functions:

(1) Solve the first-order cell functions $N_{{\alpha }_1}\left(\mathbf{y}\right)$ and $N_{{\alpha }_{1}km}\left(\mathbf{y}\right)$ based on Eqs. (14) and (15).

(2) Calculate the homogenized coefficients $G_{ij}^0$, $C_{ijkl}^0$, and $q_{\mathrm{M}}^0$ from Eq (17).

(3) Solve the second-order cell functions $N_{{\alpha }_1{\alpha }_2}\left(\mathbf{y}\right)$, $M_{{\alpha }_1{\alpha }_2}\left(\mathbf{y}\right)$, $N_{{\alpha }_1{\alpha }_{2}km}\left(\mathbf{y}\right)$, and $L_{{\alpha }_{1}k}\left(\mathbf{y}\right)$.

3. Computations of the homogenized solutions $w_0$ and ${\beta }_i^0$.

4. Assembly of the first-order and second-order multiscale approximations $w^{\epsilon ,1}$, $w^{\epsilon ,2}$, ${\beta }_i^{\epsilon ,1}$, and ${\beta }_i^{\epsilon ,2}$ from Eq (24) in an incremental way.

5. Postprocessing: Compute the first-order and second-order multiscale approximate shear stresses: $V_i^0$, $V_i^{\epsilon ,1}$, and $V_i^{\epsilon ,2}$ based on Eq (25).

As is well known, the effects of shear deformations become negligible as the thickness decreases. To address this, consistent interpolation element (CIE) is employed to solve the plate problem.

To evaluate the accuracy of the proposed asymptotic model, the relative errors are defined by

$$e_w^{\epsilon ,i}={\displaystyle \frac{{∥w^{\epsilon }-w^{\epsilon ,i}∥}_{L^2\left(\mathrm{\Omega }\right)}}{{∥w^{\epsilon }∥}_{L^2\left(\mathrm{\Omega }\right)}}},e_{{\beta }_1}^{\epsilon ,i}={\displaystyle \frac{{∥{\beta }_1^{\epsilon }-{\beta }_1^{\epsilon ,i}∥}_{L^2\left(\mathrm{\Omega }\right)}}{{∥{\beta }_1^{\epsilon }∥}_{L^2\left(\mathrm{\Omega }\right)}}},e_{{\beta }_2}^{\epsilon ,i}={\displaystyle \frac{{∥{\beta }_2^{\epsilon }-{\beta }_2^{\epsilon ,i}∥}_{L^2\left(\mathrm{\Omega }\right)}}{{∥{\beta }_2^{\epsilon }∥}_{L^2\left(\mathrm{\Omega }\right)}}},i=0,1,2,$$

(26)

where $w^{\epsilon ,0}=w_0$, ${\beta }_1^{\epsilon ,0}={\beta }_1^0$, and ${\beta }_2^{\epsilon ,0}={\beta }_2^0$.

5. Numerical Examples

5.1. Multi-Scale Approximation for Plate Problem

The SOTS asymptotic expansion method is applied to Mindlin plates made of various composite materials. Several numerical examples are provided to test the accuracy of the proposed algorithm.

5.1.1. Composite Mindlin Plate

Consider a square plate of length $L=1\mathrm{m}$ with thickness $h=0.2$m, clamped on all sides and subjected to a uniform load $q_{\mathrm{M}}=1.0\times {10}^6$ kN. The macro region and representative cell are shown in Figure 3 (a) and (b), respectively. Let material $a_1$ be carbon steel with $E^{\epsilon }=200\mathrm{GPa}$ and ${\nu }^{\epsilon }=0.25$, and material $a_2$ be low-density polyethylene(LDPE) with $E^{\epsilon }=0.3\mathrm{GPa}$ and ${\nu }^{\epsilon }=0.49$. The shear correction factor $k_s=0.833$ is adopted for both materials. $\epsilon$ is chosen as $2^{-k}$, $k=2,..,7$. The mesh information for the classical FE computations and SOTS approximation is listed in Table 1. It is evident that the classical computations take a significant amount of computational time while the proposed SOTS asymptotic method requires calculations only on a single cell and a homogenization domain with much less CPU time and also memory requirement. The relative $L^2$-errors of the approximation solutions are displayed in Table 2, which can be seen that the errors decreases as $\epsilon$ becomes small.

The convergence curves of deflection and rotations with respect to $\epsilon$ are presented in Figure 4. It is observed that the convergence orders for $w^{\epsilon }$ and ${\beta }_i^{\epsilon }$ are approximately $\mathcal{O}\left({\epsilon }^2\right)$ and $\mathcal{O}\left({\epsilon }^1\right)$, respectively. Due to all the cell functions and homogenized solutions have to be mapped into the fine meshes with different $\epsilon$. The interpolations errors between the cell mesh, homogenized mesh and the fine meshes as well as the numerical differential errors for computing the high-order derivatives of the homogenized solutions influence the relative errors computations, which to some extent conceal the convergence order of the asymptotic models.

Figure 5 and Figure 6 show the computed deflection $w^{\epsilon }$ and rotations ${\beta }_i^{\epsilon }$ when $\epsilon =1/8$. It is indicated that adding the second-order approximations $w_2$ provides a better approximation of $w^{\epsilon }$ due to the pronounced changes in deflection w. The second-order correctors ${\beta }_1^2$ and ${\beta }_2^2$ show a slight improvement over the first-order correctors ${\beta }_1^1$ and ${\beta }_2^1$, because the changes in rotation ${\beta }_1^{\epsilon }$ and ${\beta }_2^{\epsilon }$ are less obvious. Overall, the SOTS asymptotic expansions offer a more accurate representation of oscillatory behavior compared to the first-order solutions.

The approximation of the shear stresses can be obtained when $\epsilon ={\displaystyle \frac{1}{8}}$, which is shown in Figure 7. It can be observed that the oscillations are more remarkable for the stress distributions than the displacement field (Figure 7 (a) and (e)), the homogenized solutions in Figure 7 (b) and (f) only give the macroscopic behavior and the first-order solutions in Figure 7 (c) and (g) are not capable of capturing the peak of the stress within the constituent material LDPE. Therefore, the second-order correctors are necessary to be defined and appended in the asymptotic expansions to give the satisfactory distributions, which are shown in Figure 7 (d) and (h).

Figure 8 illustrates the approximations of shear stresses along the diagonal line, and it is evident that incorporating second-order correctors provides a more accurate characterization of oscillations compared to both homogenized solutions and first-order approximations.

5.1.2. Convergence for Kirchhoff Plate

Next, we check if the asymptotic model is also effective to describe the behavior of the Kirchhoff plate. Consider also a clamped square plate with a uniform load $q_{\mathrm{M}}=1.0\times {10}^3$ kN. The material parameter in $a_1$ and $a_2$ is the same as in [24]. For the macroscopic domain and each microscopic cell, divide the domain into 80 × 80 squares and 50 triangles for calculations.

To quantify the convergence behavior of the composite plate, $w_K$, ${\beta }_{1,K}$, and ${\beta }_{2,K}$ are computed using the discrete Kirchhoff theory. The computed relative $L^2$-errors are given by:

$$e_{w,K}^{\epsilon ,i}={\displaystyle \frac{{∥w^{\epsilon ,i}-w_K∥}_{L^2\left(\mathrm{\Omega }\right)}}{{∥w_K∥}_{L^2\left(\mathrm{\Omega }\right)}}},e_{{\beta }_1,K}^{\epsilon ,i}={\displaystyle \frac{{∥{\beta }_1^{\epsilon ,i}-{\beta }_{1,K}∥}_{L^2\left(\mathrm{\Omega }\right)}}{{∥{\beta }_{1,K}∥}_{L^2\left(\mathrm{\Omega }\right)}}},e_{{\beta }_2,K}^{\epsilon ,i}={\displaystyle \frac{{∥{\beta }_2^{\epsilon ,i}-{\beta }_{2,K}∥}_{L^2\left(\mathrm{\Omega }\right)}}{{∥{\beta }_{2,K}∥}_{L^2\left(\mathrm{\Omega }\right)}}}.$$

(27)

Let $h/L$ be 0.2, 0.1, and 0.035, respectively. Perform traditional FE computation for the Kirchhoff plate and compare these with asymptotic solutions. Although the plate model exhibits no significant shear locking, implementation of the CIE method is still required to prevent potential locking effects. Table 3, Table 4, and Table 5 present the relative $L^2$-errors for different $\epsilon$. It is observed that ${\beta }_1^{\epsilon }$ and ${\beta }_2^{\epsilon }$ of the Mindlin plate converge to ${\beta }_{1,K}$ and ${\beta }_{2,K}$, while $w^{\epsilon }$ approaches $w_K$ as h decreases. The deflection $w_K$ can be approximated by $w_0$ when h is not small, whereas $w^{\epsilon ,2}$ provides a more accurate approximation when h decreases. Due to the relatively smooth variation in $w_K$, $w^{\epsilon ,2}$ shows only a limited improvement in accuracy compared with $w_0$ and $w^{\epsilon ,1}$. ${\beta }_{1,K}$ and ${\beta }_{2,K}$ are best approximated by appending ${\beta }_1^2$ and ${\beta }_2^2$, which can also be shown in Figure 9 and Figure 10.

A comparison with the numerical example in [24] is performed to validate the rationality of the approach proposed in this paper. The unit cell has a height of 1 cm and in-plane dimensions $a_2$ of Y, measuring 0.3 × 0.3 cm², with a square inclusion at its center. Under the same conditions as in Example 1: Free Vibration of a CCCC Single-Layer Periodic Plate, the homogenized stiffness matrix $D^0$ can be

$$D^0=\left[\begin{array}{ccc}7.5074& 2.0974& 0\\ 2.0974& 7.5074& 0\\ 0& 0& 2.7050\end{array}\right]\times {10}^3\mathrm{N}·\mathrm{m},$$

(28)

which is closely matches the homogenized stiffness matrix $D^H$ given in [24] as follows:

$$D^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}·\mathrm{m}.$$

(29)

5.2. Multi-Scale Approximation for Beam Problem

The SOTS method is also utilized to analyze the composite Timoshenko beam.

5.2.1. Composite Timoshenko Beam

Consider a composite beam of unit length with width $b=0.1\mathrm{m}$ and height $h=0.2\mathrm{m}$, which is clamped at both ends and subjected to a uniform load $q_{\mathrm{T}}=1.0\times {10}^5\mathrm{kN}$. The macro and cell domains are shown in Figure 11 (a) and (b), where ${\mathrm{\Omega }}_a$(${\mathrm{\Omega }}_b$) is occupied by the material $y_a$($y_b$). Choose $E^{\epsilon }=17\mathrm{GPa}$ and ${\nu }^{\epsilon }=0.42$ for the plumbum ($y_a$) in ${\mathrm{\Omega }}_a$, and $E^{\epsilon }=206\mathrm{GPa}$ and ${\nu }^{\epsilon }=0.3$ for the alloy steel ($y_b$) in ${\mathrm{\Omega }}_b$. The macroscopic domain $\mathrm{\Omega }$ is discretized into 100 finite elements for computation. Let $\epsilon$ be $2^{-k}$, $k=2,...,8$. The computed relative $L^2$-errors are shown in Table 6.

It is shown that the SOTS computation also offers improved precision for approximating both $w^{\epsilon }$ and ${\beta }^{\epsilon }$ with less computational effort compared to the general FE method.

The asymptotic approximations for $\epsilon ={\displaystyle \frac{1}{2}}$, ${\displaystyle \frac{1}{4}}$, and ${\displaystyle \frac{1}{8}}$ are presented in Figure 12, with local detail graphs within. Compared to $\epsilon ={\displaystyle \frac{1}{2}}$, $w^{\epsilon ,2}$ and ${\beta }^{\epsilon ,2}$ are gradually closer to $w^{\epsilon }$ and ${\beta }^{\epsilon }$ for $\epsilon ={\displaystyle \frac{1}{4}}$ and ${\displaystyle \frac{1}{8}}$. As shown in the figures, the approximation for $\epsilon ={\displaystyle \frac{1}{8}}$ is more effective compared to $\epsilon ={\displaystyle \frac{1}{4}}$ with more details captured.

The converging behaviors of the deflection and the rotation with respect to $\epsilon$ are shown in Figure 13. It can be observed from the figure that the convergence rates for the $w^{\epsilon }$ and ${\beta }^{\epsilon }$ are approximately $\mathcal{O}\left({\epsilon }^2\right)$, whereas the convergence order of the homogenized result is approximately $\mathcal{O}\left({\epsilon }^1\right)$.

5.2.2. Convergence for Euler Beam

Choose $q_{\mathrm{T}}=1.0\times {10}^3$$\mathrm{kN}$. As the thickness h decreases, the beam model requires finer mesh discretization. In the macroscopic domain, the domain is divided into 400 mesh elements for calculations, and other parameters are selected similarly to those in Section 5.2.1. Define $w_E$ and ${\beta }_E$ as the deflection and rotation of the Euler beam obtained by the general FE method. The relative $L^2$-errors $e_{w,E}^{\epsilon ,i}$ and $e_{\beta ,E}^{\epsilon ,i}$ are defined similarly to (27).

Compare the solutions of the numerical method with the solutions $w_E$ and ${\beta }_E$ of the Euler beam using the general finite element method. Since the rotation ${\beta }^{\epsilon }$ of the Timoshenko beam remains close to the rotation ${\beta }_E$ for $h=0.1\mathrm{m}$, $0.035\mathrm{m}$, and $0.02\mathrm{m}$, while the deflection $w^{\epsilon }$ of the Timoshenko beam approximates $w_E$ more accurately when h is small, as shown in Figure 14 (a)-(f). It is shown that approximating the deflection $w_E$ by appending the second-order corrector $w_2$ is more effective than using the first-order corrector $w_1$ and the homogenized solution $w_0$, when h is small as $h=0.035\mathrm{m}$ and $0.02\mathrm{m}$. The microscopic homogenized coefficients used to obtain the homogenized Euler beam are displayed as follows:

$$E^0=24503,{\nu }^0=0.415236,G^0=8,656.8601.$$

(30)

The relative $L^2$-errors of the homogenized Euler beam when approximating the composite Euler beam are displayed in Table 7. As $\epsilon$ gradually decreases, the relative errors are reduced, and approximations of both deflection and rotation improve.

Additionally, the rotation ${\beta }_E$ is better approximated by adding the second-order corrector ${\beta }_2$ than using the homogenized solution ${\beta }_0$ or appending the first-order corrector ${\beta }_1$ for different values of h. This is detailed in Table 8, Table 9, and Table 10. As h decreases from 0.1 m to 0.02 m, the composite Timoshenko beam gradually approximates the composite Euler beam more closely.

When the thickness h reaches a minimum of 0.02 m, no significant shear locking phenomenon is observed. However, when h is less than 0.001 m, pronounced shear locking occurs, necessitating the application of RIE method to address this numerical issue.

6. Conclusions

A two-scale asymptotic framework is established for composite Mindlin plates and Timoshenko beams, demonstrating that the second-order correctors are crucial for accurately capturing stress oscillations at material interfaces. When h decreases, the asymptotic model developed can effectively approximate the Kirchhoff plate and Euler beam with composite configurations, respectively, providing an alternative way to simulate complicated thin plate and beam structures. The RIE and CIE are employed in the finite element computations of the homogenization models to mitigate shear locking, which can give an accurate description of the displacement and stress distributions of the plates and beams.

Most importantly, the proposed method offers significant potential for future work. It can be directly extended to model sustainable composites with valorised manufacturability. This establishes our asymptotic framework as a general and computationally efficient approach for predicting the multi-physical response of next-generation composite plates and beams, thereby contributing to the development of tailored material systems aligned with the goals of modern industrial applications.