A novelty solution for orthotropic composite plate based on physics informed neural network

Abstract

The modeling of cross-ply composite laminates using numerical methods has been a difficult task, leading to the development of various finite element method and other analytical solutions. However, as materials science advances, this problem has become more complex, presenting new challenges that require reliable and novel approaches. In this study, we propose the utilization of machine learning, specifically physics informed neural networks (PINN), for the first time to examine the behavior of composite plate. By solving a system of partial differential equations derived from the virtual work equilibrium principle, PINN are employed as a method to solve these equations using a generalized strong-form approach. To address the issue of imbalanced loss functions, we also propose adjusting the loss function in this research. Once trained, PINN serve as a surrogate model capable of predicting displacements and stresses in cross-ply composite laminates. To demonstrate the effectiveness and reliability of PINN, we investigate two examples of laminates with different material distributions and boundary conditions including boundary conditions on displacement and boundary conditions on stress, comparing the results with the benchmark Navier solution. The research and obtained results showcase the performance and accuracy of PINN, highlighting their potential as a surrogate model for solving problems related to cross-ply composite laminates.

Graphical abstract

Keywords

PINN / loss function / laminate plates / Navier solution

1 Introduction

Because of their remarkable mechanical characteristics and lightweight design, laminated composite materials have drawn a lot of interest from engineers for a variety of technical applications [1–4]. The study examined the combination of flow dynamics and data-driven modeling utilizing physics informed neural networks (PINN). By merging empirical data with theoretical frameworks, PINN successfully tackled inverse problems across multiple fields, such as three-dimensional wake phenomena, high-speed aerodynamic flows, and biomedical fluid dynamics [5].

Analyzing the static behavior of laminated plates is crucial for understanding their structural performance and optimizing their design. There have been numerous studies developed to enrich computational methods, including various numerical solutions such as finite element method (FEM) [6–8]. However, it can be difficult to adequately capture the complicated behavior of laminated plates using classic analytical and numerical approaches, like finite element analysis, especially when dealing with complex boundary conditions and material distributions [9,10]. Hence, research endeavors to explore alternative and effective methods, instead of relying solely on traditional approaches, continue to be pursued and developed. In the past decade, the rapid progress of artificial neural networks (NNs) and optimization techniques has facilitated the successful resolution of various engineering challenges, yielding impressive accuracy and promising results [11–13]. The remarkable predictive capability of artificial NNs, as highlighted by the universal approximation theorem, is the key factor behind their effectiveness [14]. With the proliferation of data collected from diverse sources and the breakthroughs achieved in data-driven deep learning, researchers have been inspired to extend its robustness to the realm of applied mechanics. This extension has led to numerous successful applications in the field of structural engineering [15–18]. Raissi et al. [19] introduced an advanced deep learning approach derived from the original PINN to solve partial differential equations (PDEs). PINN combine the power of NNs with the principles of physics to provide accurate and efficient solutions. By incorporating physical laws and data into the learning process, PINN can effectively capture the underlying physics of the problem. Most engineering problems in contemporary times are built upon the foundations of mathematics and physics. As a result, the equations that can describe the behaviors of these problems often stem from differential equations. As a result, the solution of differential equations plays a critical role in addressing practical problems [20,21]. Utilizing PINN for solving PDEs offers multiple advantages: physical laws and data can be incorporated [22]; complex problems can be handled with flexibility [23]; handling complex geometries [24,25]; and parallelization [26]. Because of these benefits, PINN are a viable method for solving PDEs accurately and quickly in a variety of scientific and technical fields [27–32]. In recent years, the use of PINN to solve mechanical problems has gained significant attention [5,30,33]. PINN combine the power of NNs with the principles of physics to provide accurate and efficient solutions. With PINN, it is possible to build models based on physical principles and governing equations to simulate and predict mechanical phenomena in systems. PINN can learn from available data, including boundary conditions and known solutions, to find accurate functional representations of the spatial and temporal variables of the mechanical quantities [34–37]. PINN have been used as an alternative method to solve various engineering problems.

Machine learning methods have surfaced as viable substitutes in recent years for resolving challenging engineering issues. The application of PINN to the static analysis of laminated plates has demonstrated significant promise among these methods. PINN offer precise and effective solutions by fusing the capabilities of NNs with the laws of physics. PINN are highly suited for static analysis of laminated plates because they can effectively capture the underlying physics of the problem by adding physical principles and governing equations into the learning process, several studies have been conducted, and promising results have been achieved, as PINN can accurately predict the plate’s displacements and stresses with a high level of accuracy. The initial studies on the use of PINN for Kirchhoff plate problems were introduced by Ref. [38]. In this research, a deep collocation method (DCM) for thin plate bending problems was presented, leveraging computational graphs and backpropagation algorithms from deep learning. This method, based on a feedforward of deep neural network (DNN), addressed the challenges of C1 continuity in traditional mesh-based methods by approximating continuous transversal deflection. The approach effectively minimized the governing PDE and boundary conditions, proving suitable for the bending analysis of Kirchhoff plates with various geometries. Zhuang et al. [39] introduced a deep autoencoder based energy method (DAEM) for analyzing the bending, vibration, and buckling behavior of Kirchhoff plates. This approach capitalized on the model’s higher-order continuity and seamlessly integrated a deep autoencoder following the principle of minimizing total potential energy into a unified framework, forming an unsupervised feature learning technique. As a specialized form of a feedforward of DNN, the DAEM served as a function approximator. Leveraging its powerful feature extraction capabilities, the model efficiently recognized intricate patterns throughout the entire energy system. This included analyzing field variables, natural frequencies, and critical buckling load factors, all of which were central to the study. The primary goal of the optimization process was to achieve a state where the total potential energy reached its minimum, aligning with fundamental principles in structural analysis.

Li et al. [22] evaluated the robustness of PINN by analyzing four unique loading conditions: non-uniformly distributed tensile forces in in-plane stretching, central-hole tension under in-plane loading, out-of-plane displacement, and compression-driven buckling. Three distinct loss function formulations were examined: 1) entirely data-driven, 2) governed by PDEs, and 3) energy-based. By benchmarking the results against finite element simulations, it was found that all three methods effectively captured the elastic deformation of plates with a commendable degree of accuracy. Bastek and Kochmann [40] applied PINN to estimate the small-strain behavior of arbitrarily curved shells. In this research, the widely recognized Scordelis-Lo roof model was used to evaluate the effectiveness of PINN in comparison with the FEM. The findings demonstrate that when the equations are expressed in their weak formulation, the PINN can precisely determine the solution field across all three benchmark cases. Samaniego et al. [41] introduced an energy-based method for solving PDEs in computational mechanics through machine learning. In this study, DNNs were effectively utilized in computational mechanics to establish an approximation space for addressing various important cases. The fast progression of deep learning has had a profound impact on disciplines such as solid mechanics, especially in tackling PDEs using PINN and the deep energy method (DEM) to address the limitations of DEM, which relied solely on the principle of minimum potential energy. On the basis of the principle of minimum complementary energy, Wang et al. [42] introduced the deep complementary energy method (DCEM). The computational findings revealed that DCEM surpassed DEM in terms of stress precision and computational efficiency, particularly in scenarios with intricate boundary conditions. This demonstrates the feasibility of utilizing deep learning concepts and tools to address the solution of highly relevant boundary value problems. Peng et al. [43] introduced a PINN method with a two-network strategy to solve the bending problem of thin plates with variable stiffness on an elastic foundation. Instead of directly solving the fourth-order PDE, the method transformed it into four second-order PDEs based on Kirchhoff plate theory. The numerical results were validated against FEM and literature, demonstrating the method's effectiveness, stability, and ease of applying boundary conditions.

Fallah and Aghdam [44] employed PINN to evaluate the bending response and free vibration characteristics of three-dimensional functionally graded (TDFG) beams with spatially varying material properties. The governing motion equations are addressed through the PINN approach, and the outcomes are benchmarked against reference solutions. Additionally, the research examines the influence of material gradation, elastic foundation support, porosity levels, and different porosity distribution patterns on the bending performance and natural frequencies of TDFG beams. Wang and Thai [45] developed a PINN framework for predicting the bending behavior of clamped laminated composite plates. Comparisons with analytical and finite element results show that the PINN framework achieves accurate predictions with improved computational efficiency. Vahab et al. [46] used PINN to combine with Airy stress functions and Fourier series to optimize solutions for challenging biharmonic problems in elasticity and elastic plate theory. By integrating classical analytical methods, efficient NNs with minimal parameters are constructed, offering highly accurate and fast evaluations. Incorporating Airy stress functions in the feature space notably improves the accuracy of PINN solutions for biharmonic problems. Mishra et al. [47] investigated the stability loss of a laminated composite skew plate utilizing a C₀ finite element framework grounded in advanced shear deformation theory. Furtado et al. [48] utilized data-driven techniques to approximate the statistical design allowable of composite laminates. Four distinct machine learning approaches including artificial NNs, XGBoost, Random Forests and Gaussian Processes were employed to forecast notched strength and its distribution while incorporating variability in material properties and geometric attributes. These models demonstrated exceptional precision, accurately representing the design space with minimal overfitting. The results provide important perspectives on predicting failure behaviors and ultimate strength, significantly cutting down computational costs in establishing design allowable for composite laminates. A recent investigation presented Kolmogorov-Arnold-Informed Neural Networks (KINNs) [49] as an advancement over conventional multilayer perceptron-based PINN for addressing solid mechanics challenges. KINNs enhanced accuracy and accelerated convergence, especially when solving PDEs in computational solid mechanics, though their effectiveness was limited for intricate geometries. This methodology demonstrated potential for delivering more efficient and precise solutions within the domain.

This study is the first to tackle this intricate issue. The paper outlines the problem formulation, the training process of the PINN model, and its validation through comparisons with analytical and experimental data. Numerical case studies will illustrate the efficiency and practical applications of the proposed method, highlighting the benefits of employing PINN for the static analysis of laminated plates

2 The classical laminated plate theory

2.1 Strain–displacement relation

The CLPT builds upon the classical plate theory to accommodate composite laminates [50]. Within CLPT, Kirchhoff’s assumptions are adopted as follows. 1) Straight lines perpendicular to the mid-surface of the plate remain perpendicular after deformation and are therefore referred to as transverse normal. 2) The transverse normal stay straight even after undergoing deformation. 3) The plate’s thickness remains unchanged post-deformation.

Consider a plate laminate has a thickness h and are made up of Nanisotropic layers with material orientations \left(x_1^k,x_2^k,x_3^k\right) of the kth lamina secured by an angle θth to the system of laminate coordinate. Assume that z axis is positive downward given in Fig.1. The location of kth lamina is determined based on two points \left({\mathit{\text{z}}}_k,{\mathit{\text{z}}}_{k+1}\right) as shown in Fig.1.

Based on the Kirchhoff assumptions, the displacements \left({λ̶¯}_{\mathrm{u}},{λ̶¯}_{\mathrm{v}},{λ̶¯}_{\mathrm{w}}\right) at any points \left(x,y\right) in laminate plate using CLPT as shown in Eq. (1).

(1)\begin{array}{rl}& {λ̶¯}_{\mathrm{u}}={λ̶˜}_{{\mathrm{u}}_0}-\mathit{\text{z}}\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }x},\\ & {λ̶¯}_{\mathrm{v}}={λ̶˜}_{{\mathrm{v}}_0}-\mathit{\text{z}}\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }y},\\ & {λ̶¯}_{\mathrm{w}}={λ̶˜}_{{\mathrm{w}}_0},\end{array}

where ({λ̶˜}_{{\mathrm{u}}_0},{λ̶˜}_{{\mathrm{v}}_0},{λ̶˜}_{{\mathrm{w}}_0}) are the corresponding displacement at the mid-plane of laminated plate at coordinate \left(x,y,\mathit{\text{z}}\right) of displacement {λ̶¯}_{\mathrm{u}}, {λ̶¯}_{\mathrm{v}} and {λ̶¯}_{\mathrm{w}}. When we determine the values of ( {λ̶˜}_{{\mathrm{u}}_0},{λ̶˜}_{{\mathrm{v}}_0},{λ̶˜}_{{\mathrm{w}}_0}), we will specify the displacement at point having coordinate \left(x,y\right) in the domain of laminated plate. Because the second-order derivatives are negligible in the deformation equation, they are thus ignored. The deformation in the global coordinate \left(x,y,\mathit{\text{z}}\right) is given by the following Eq. (2).

(2)\begin{array}{rl}& {\epsilon }_{xx}=\frac{\mathrm{\partial }\left({λ̶¯}_{\mathrm{u}}\right)}{\mathrm{\partial }x}+\frac{1}{2}{\left(\frac{\mathrm{\partial }\left({λ̶¯}_{\mathrm{w}}\right)}{\mathrm{\partial }x}\right)}^2,{\epsilon }_{yy}=\frac{\mathrm{\partial }\left({λ̶¯}_{\mathrm{v}}\right)}{\mathrm{\partial }y}+\frac{1}{2}{\left(\frac{\mathrm{\partial }\left({λ̶¯}_{\mathrm{w}}\right)}{\mathrm{\partial }y}\right)}^2,\\ & {\epsilon }_{xy}=\frac{1}{2}\left(\frac{\mathrm{\partial }\left({λ̶¯}_{\mathrm{u}}\right)}{\mathrm{\partial }y}+\frac{\mathrm{\partial }\left({λ̶¯}_{\mathrm{v}}\right)}{\mathrm{\partial }x}+\frac{\mathrm{\partial }\left({λ̶¯}_{\mathrm{w}}\right)}{\mathrm{\partial }x}\frac{\mathrm{\partial }\left({λ̶¯}_{\mathrm{w}}\right)}{\mathrm{\partial }y}\right),\\ & {\epsilon }_{x\mathit{\text{z}}}=\frac{1}{2}\left(\frac{\mathrm{\partial }\left({λ̶¯}_{\mathrm{u}}\right)}{\mathrm{\partial }\mathit{\text{z}}}+\frac{\mathrm{\partial }\left({λ̶¯}_{\mathrm{w}}\right)}{\mathrm{\partial }\mathit{\text{z}}}\right),{\epsilon }_{y\mathit{\text{z}}}=\frac{1}{2}\left(\frac{\mathrm{\partial }\left({λ̶¯}_{\mathrm{v}}\right)}{\mathrm{\partial }\mathit{\text{z}}}+\frac{\mathrm{\partial }\left({λ̶¯}_{\mathrm{w}}\right)}{\mathrm{\partial }y}\right),\\ & {\epsilon }_{\mathit{\text{zz}}}=\frac{\mathrm{\partial }\left({λ̶¯}_{\mathrm{w}}\right)}{\mathrm{\partial }\mathit{\text{z}}}.\end{array}

For the assumed displacement field, \mathrm{\partial }\left({λ̶¯}_{\mathrm{w}}\right)/\mathrm{\partial }\mathit{\text{z}}=0. Equation (2) can be rewritten by substituting Eq. (1) into Eq. (2). We will obtain the deformation equation related to the values \left({λ̶˜}_{{\mathrm{u}}_0},{λ̶˜}_{{\mathrm{v}}_0},{λ̶˜}_{{\mathrm{w}}_0}\right)as presented in Eq. (3).

(3)\begin{array}{rl}& {\epsilon }_{xx}=-\mathit{\text{z}}\frac{{\mathrm{\partial }}^2\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }{x}^2}+\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{u}}_0}\right)}{\mathrm{\partial }x}+\frac{1}{2}{\left(\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }x}\right)}^2,\\ & {\epsilon }_{xy}=-\mathit{\text{z}}\frac{{\mathrm{\partial }}^2\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }x\mathrm{\partial }y}+\frac{1}{2}\left(\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{u}}_0}\right)}{\mathrm{\partial }y}+\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{v}}_0}\right)}{\mathrm{\partial }x}+\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }x}\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }y}\right),\\ & {\epsilon }_{yy}=-\mathit{\text{z}}\frac{{\mathrm{\partial }}^2\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }{y}^2}+\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{v}}_0}\right)}{\mathrm{\partial }y}+\frac{1}{2}{\left(\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }y}\right)}^2,\\ & {\epsilon }_{x\mathit{\text{z}}}=\frac{1}{2}\left(-\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }x}+\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }x}\right)=0,\\ & {\epsilon }_{y\mathit{\text{z}}}=\frac{1}{2}\left(-\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }y}+\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }y}\right)=0,\\ & {\epsilon }_{\mathit{\text{zz}}}=0.\end{array}

Equation (3) can be expressed concisely as Eq. (4)

(4)\left\{\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\Gamma }_{xy}\end{array}\right\}=\left\{\begin{array}{c}{\epsilon }_{xx}^{\left(0\right)}\\ {\epsilon }_{yy}^{\left(0\right)}\\ {\Gamma }_{yy}^{\left(0\right)}\end{array}\right\}+\mathit{\text{z}}\left\{\begin{array}{c}{\epsilon }_{xx}^{\left(1\right)}\\ {\epsilon }_{yy}^{\left(1\right)}\\ {\Gamma }_{yy}^{\left(1\right)}\end{array}\right\},

where \left({\epsilon }_{xx}^{\left(0\right)},{\epsilon }_{yy}^{\left(0\right)},{\Gamma }_{xy}^{\left(0\right)}\right) are the membrane strains and the other one \left({\epsilon }_{xx}^{\left(1\right)},{\epsilon }_{yy}^{\left(1\right)},{\Gamma }_{xy}^{\left(1\right)}\right) are the flexural bending strains. They are shown in explicit form as Eq. (5).

(5)\begin{array}{rl}& \left\{\begin{array}{c}{\epsilon }_{xx}^{\left(0\right)}\\ {\epsilon }_{yy}^{\left(0\right)}\\ {\Gamma }_{yy}^{\left(0\right)}\end{array}\right\}=\left\{\begin{array}{c}{\displaystyle \frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{u}}_0}\right)}{\mathrm{\partial }x}}+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }x}}\right)}^2\\ {\displaystyle \frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{v}}_0}\right)}{\mathrm{\partial }y}}+{\displaystyle \frac{1}{2}}{\left({\displaystyle \frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }y}}\right)}^2\\ {\displaystyle \frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }x}}{\displaystyle \frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }y}}+{\displaystyle \frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{u}}_0}\right)}{\mathrm{\partial }y}}+{\displaystyle \frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{v}}_0}\right)}{\mathrm{\partial }x}}\end{array}\right\},\\ & \left\{\begin{array}{c}{\epsilon }_{xx}^{\left(1\right)}\\ {\epsilon }_{yy}^{\left(1\right)}\\ {\Gamma }_{yy}^{\left(1\right)}\end{array}\right\}=\left\{\begin{array}{c}-{\displaystyle \frac{{\mathrm{\partial }}^2\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }{x}^2}}\\ -{\displaystyle \frac{{\mathrm{\partial }}^2\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }{y}^2}}\\ -2{\displaystyle \frac{{\mathrm{\partial }}^2\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }x\mathrm{\partial }y}}\end{array}\right\}.\end{array}

2.2 Lamina stress–strain relation

In the CLPT, the strain components ({\epsilon }_{\mathit{\text{zz}}},{\epsilon }_{x\mathit{\text{z}}} ,{\epsilon }_{y\mathit{\text{z}}} ) are assumed to be zero. As a result, the shear stresses \left({\sigma }_{x\mathit{\text{z}}},{\sigma }_{y\mathit{\text{z}}}\right) are also to be zero. To establish the governing equations of a laminate in the global coordinate related to the material coordinates \left(x_1,x_2,x_3\right), We assume that \phi is the angle between the x_1 axis of the material coordinate system and the xaxis of the global coordinate system, as shown in Fig.2. The relationship between the stresses in \left(x,y,\mathit{\text{z}}\right) and the material coordinates \left(x_1,x_2,x_3\right) of the k^th orthotropic lamina are shown in Eq. (6).

(6){\left\{\begin{array}{c}{\sigma }_1\\ {\sigma }_2\\ {\sigma }_6\end{array}\right\}}^{\left(k\right)}={\left[\begin{array}{ccc}{R}_{11}& R_{12}& 0\\ R_{12}& R_{22}& 0\\ 0& 0& R_{66}\end{array}\right]}^{\left(k\right)}\left\{\begin{array}{c}{\epsilon }_1\\ {\epsilon }_2\\ {\epsilon }_6\end{array}\right\},

where \left({\sigma }_1={\sigma }_{11},{\sigma }_2={\sigma }_{22},{\sigma }_6={\sigma }_{12}\right), and similarly for the deformation components. The symbol R_{ij}^{\left(k\right)} are the plane stress-reduced stiffness. They are engineering constants registered at kth layer given in Eq. (7).

(7)\begin{array}{rl}& R_{11}=\frac{E_1}{1-v_{12}{v}_{21}},R_{12}=\frac{v_{12}{E}_2}{1-v_{12}{v}_{21}}=\frac{v_{21}{E}_1}{1-v_{12}{v}_{21}},\\ & R_{22}=\frac{E_2}{1-v_{12}{v}_{21}},R_{22}=G_{12},\end{array}

where \left(E_1,E_2,G_{12},v_{12},v_{21}\right) are the material parameters. Because the laminate is created by many layers, therefore, the constitutive relations for each layer need to be converted into the global coordinate \left(x,y,\mathit{\text{z}}\right).The relationship between the stresses \left({\sigma }_{xx},{\sigma }_{yy},{\sigma }_{xy}\right) and strain \left({\epsilon }_{xx},{\epsilon }_{yy},{\epsilon }_{xy}\right) as shown in Eq. (8).

(8){\left\{\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\sigma }_{xy}\end{array}\right\}}^{\left(k\right)}={\left[\begin{array}{ccc}{\overline{R}}_{11}& {\overline{R}}_{12}& {\overline{R}}_{16}\\ {\overline{R}}_{12}& {\overline{R}}_{22}& {\overline{R}}_{26}\\ {\overline{R}}_{16}& {\overline{R}}_{26}& {\overline{R}}_{66}\end{array}\right]}^{\left(k\right)}\left\{\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\epsilon }_{xy}\end{array}\right\},

where {\overline{R}}_{11}^{\left(k\right)}, {\overline{R}}_{12}^{\left(k\right)}, {\overline{R}}_{22}^{\left(k\right)}, {\overline{R}}_{66}^{\left(k\right)} are the plan stress-reduced stiffness of the kth lamina referred to the material coordinate system oriented by an angle \phi and can be expressed as Eq. (9).

(9)\begin{array}{rl}& {\overline{R}}_{11}^k=R_{11}^k{\mathrm{c}\mathrm{o}\mathrm{s}}^4\phi +2(R_{12}^k+2Q_{66}^k){\mathrm{s}\mathrm{i}\mathrm{n}}^2\phi {\mathrm{c}\mathrm{o}\mathrm{s}}^2\phi +R_{22}^k{\mathrm{s}\mathrm{i}\mathrm{n}}^4\phi ,\\ & {\overline{R}}_{12}^k=(R_{11}^k+R_{22}^k-4R_{66}^k){\mathrm{s}\mathrm{i}\mathrm{n}}^2\phi {\mathrm{c}\mathrm{o}\mathrm{s}}^2\phi +R_{12}^k({\mathrm{s}\mathrm{i}\mathrm{n}}^4\phi +{\mathrm{c}\mathrm{o}\mathrm{s}}^4\phi ),\\ & {\overline{R}}_{22}^k=R_{11}^k{\mathrm{s}\mathrm{i}\mathrm{n}}^4\phi +2(R_{12}^k+2R_{66}^k){\mathrm{s}\mathrm{i}\mathrm{n}}^2\phi {\mathrm{c}\mathrm{o}\mathrm{s}}^2\phi +R_{22}^k{\mathrm{c}\mathrm{o}\mathrm{s}}^4\phi ,\\ & {\overline{R}}_{16}^k=(R_{11}^k-R_{12}^k-2{\overline{R}}_{66}^k)\mathrm{s}\mathrm{i}\mathrm{n}\phi {\mathrm{c}\mathrm{o}\mathrm{s}}^3\phi +(R_{12}^k-R_{22}^k+2R_{66}^k){\mathrm{s}\mathrm{i}\mathrm{n}}^3\phi \mathrm{c}\mathrm{o}\mathrm{s}\phi ,\\ & {\overline{R}}_{26}^k=(R_{11}^k-R_{12}^k-2{\overline{R}}_{66}^k)\mathrm{s}\mathrm{i}\mathrm{n}\phi {\mathrm{c}\mathrm{o}\mathrm{s}}^3\phi +(Q_{12}^k-Q_{22}^k+2R_{66}^k)\mathrm{s}\mathrm{i}\mathrm{n}\phi {\mathrm{c}\mathrm{o}\mathrm{s}}^3\phi ,\\ & {\overline{R}}_{66}^k=(R_{11}^k+R_{22}^k-2R_{12}^k-2R_{66}^k){\mathrm{s}\mathrm{i}\mathrm{n}}^2\phi {\mathrm{c}\mathrm{o}\mathrm{s}}^2\phi +R_{66}^k({\mathrm{s}\mathrm{i}\mathrm{n}}^4\phi +{\mathrm{c}\mathrm{o}\mathrm{s}}^4\phi ).\end{array}

2.3 Governing equations

In this study, edge forces, thermal expansion, and electric fields are neglected, and the governing equations are formulated based on the principle of virtual displacement. Notations\left({\hat{\sigma }}_{nn},{\hat{\sigma }}_{ns},{\hat{\sigma }}_{n\mathit{\text{z}}}\right) are the stress component on the boundary {\Gamma }_{\sigma } which is a part of \Gamma. And \left(\delta {λ̶¯}_{u_{0n}},\delta {λ̶¯}_{u_{0s}}\right) are the virtual displacements according to normal and tangential directions, respectively given in Fig.3.

The equilibrium form of the virtual work principle is given in Eq. (10).

(10)0=\delta U+ \delta V,

where \delta U is the strain energy and \delta Vis the virtual work secured by an applied force. recall the Eq. (4) and the CLPT theory \delta {\epsilon }_{x\mathit{\text{z}}}=0, \delta {\epsilon }_{y\mathit{\text{z}}} =0, \delta {\epsilon }_{\mathit{\text{zz}}}=0. So, their extended forms are presented in the following Eqs. (11) and (12).

(11)\begin{array}{rl}\delta U=& \int {\mathrm{\Omega }}_0{\int }_{-h/2}^{h/2}\left({\sigma }_{xx}\delta {\epsilon }_{xx}+{\sigma }_{yy}\delta {\epsilon }_{yy}+2{\sigma }_{xy}\delta {\epsilon }_{xy}\right)\mathrm{d}x\mathrm{d}y\mathrm{d}\mathit{\text{z}}=\\ & \int {\mathrm{\Omega }}_0{\int }_{-h/2}^{h/2}\left({\sigma }_{xx}\left(\delta {\epsilon }_{xx}^{\left(0\right)}+\mathit{\text{z}}\delta {\epsilon }_{xx}^{\left(1\right)}\right)+{\sigma }_{yy}\left(\delta {\epsilon }_{yy}^{\left(0\right)}+\mathit{\text{z}}\delta {\epsilon }_{yy}^{\left(1\right)}\right)+2{\sigma }_{xy}\left({\epsilon }_{xy}^{\left(0\right)}+\mathit{\text{z}}\delta {\epsilon }_{xy}^{\left(1\right)}\right)\right)\mathrm{d}x\mathrm{d}y\mathrm{d}\mathit{\text{z}},\end{array}

(12)\begin{array}{rl}& \delta V=-\int {\mathrm{\Omega }}_0\left[q_{\mathrm{b}}\delta {λ̶¯}_{\mathrm{w}}+q_{\mathrm{t}}\delta {λ̶¯}_{\mathrm{w}}\right]\mathrm{d}x\mathrm{d}y-\int {\Gamma }_{\sigma }{\int }_{-h/2}^{h/2}\left[{\hat{\sigma }}_{nn}\delta {λ̶¯}_{u_n}+{\hat{\sigma }}_{ns}\delta {λ̶¯}_{u_s}+{\hat{\sigma }}_{n\mathit{\text{z}}}\delta {λ̶¯}_{\mathrm{w}}\right]\mathrm{d}\mathit{\text{z}}\mathrm{d}s\\ & \delta V=-\int {\mathrm{\Omega }}_0\left[q_{\mathrm{b}}+q_{\mathrm{t}}\right]\delta {λ̶˜}_{{\mathrm{w}}_0}\mathrm{d}x\mathrm{d}y-\int {\Gamma }_{\sigma }{\int }_{-h/2}^{h/2}\left[{\hat{\sigma }}_{nn}\left(\delta {λ̶˜}_{u_{0n}}-\mathit{\text{z}}\frac{\mathrm{\partial }\left(\delta {λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }n}\right)+{\hat{\sigma }}_{ns}\left(\delta {λ̶˜}_{u_{0s}}-\mathit{\text{z}}\frac{\mathrm{\partial }\left(\delta {λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }s}\right)+{\hat{\sigma }}_{n\mathit{\text{z}}}\delta {λ̶˜}_{{\mathrm{w}}_0}\right]\mathrm{d}\mathit{\text{z}}\mathrm{d}s,\end{array}

where q_{\mathrm{b}}\left(x,y\right) and q_{\mathrm{t}}\left(x,y\right) are the distributed forces at the bottom and top, respectively. {\mathrm{\Omega }}_0 denote the midplane domain of the plate. The virtual strains are determined in a manner analogous to the true strain, both of which are expressed in terms of the actual displacement field. The explicit formulation of these virtual strains can be presented in detail in Eq. (13).

(13)\begin{array}{rl}\delta {\epsilon }_{xx}^{\left(0\right)}=& \frac{\mathrm{\partial }\left(\delta {λ̶˜}_{{\mathrm{u}}_0}\right)}{\mathrm{\partial }x}+\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }x}\frac{\mathrm{\partial }\left(\delta {λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }x},\\ \delta {\epsilon }_{yy}^{\left(0\right)}=& \frac{\mathrm{\partial }\left(\delta {λ̶˜}_{{\mathrm{v}}_0}\right)}{\mathrm{\partial }y}+\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }y}\frac{\mathrm{\partial }\left(\delta {λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }y},\\ \delta {\Gamma }_{xy}^{\left(0\right)}=& \frac{\mathrm{\partial }\left(\delta {λ̶˜}_{{\mathrm{u}}_0}\right)}{\mathrm{\partial }x}+\frac{\mathrm{\partial }\left(\delta {λ̶˜}_{{\mathrm{v}}_0}\right)}{\mathrm{\partial }y}+\frac{\mathrm{\partial }\left(\delta {λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }x}\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }y}\\ & +\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }x}\frac{\mathrm{\partial }\left(\delta {λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }y},\delta {\epsilon }_{xx}^{\left(1\right)}=-\frac{{\mathrm{\partial }}^2\left(\delta {λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }{x}^2},\\ \delta {\epsilon }_{yy}^{\left(1\right)}=& -\frac{{\mathrm{\partial }}^2\left(\delta {λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }{y}^2},\delta {\Gamma }_{xy}^{\left(1\right)}=-2\frac{{\mathrm{\partial }}^2\left(\delta {λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }x\mathrm{\partial }y}.\end{array}

By replacing Eqs. (11)–(13) into Eq. (10) and integrating across the thickness of the plate, Eq. (10) is rewritten as

(14)\begin{array}{rl}0=& {\int }_{{\mathrm{\Omega }}_0}[-\left({\stackrel{x}{N}}_{xx}+{\stackrel{y}{N}}_{xy}\right)\delta {λ̶˜}_{{\mathrm{u}}_0}-\left({\stackrel{x}{N}}_{xy}+{\stackrel{y}{N}}_{yy}\right)\delta {λ̶˜}_{{\mathrm{v}}_0}-\left({\stackrel{xx}{M}}_{xx}+2{\stackrel{xy}{M}}_{xy}+{\stackrel{yy}{M}}_{yy}+\mathrm{\aleph }\left({λ̶˜}_{{\mathrm{w}}_0}\right)+q\right)\delta {λ̶˜}_{w_0}]\mathrm{d}x\mathrm{d}y\\ & +{\int }_{{\Gamma }_{\sigma }}[\left(N_{xx}{\tau }_x+N_{xy}{\tau }_y\right)\delta {λ̶˜}_{{\mathrm{u}}_0}+\left(N_{xy}{\tau }_x+N_{yy}{\tau }_y\right)\delta {λ̶˜}_{{\mathrm{v}}_0}+\left({\stackrel{x}{M}}_{xx}{\tau }_x+M_{xy}{\tau }_y+{\stackrel{x}{M}}_{yy}{\tau }_y+{\stackrel{x}{M}}_{xy}{\tau }_y+\rho \left({λ̶˜}_{{\mathrm{w}}_0}\right)\right)\delta {λ̶˜}_{{\mathrm{w}}_0}\\ & -\left(M_{xx}{\tau }_x+M_{xy}{\tau }_y\right)\frac{\mathrm{\partial }\left(\delta {λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }x}-\left(M_{xy}{\tau }_x+M_{yy}{\tau }_y\right)\frac{\mathrm{\partial }\left(\delta {λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }y}]\mathrm{d}s\\ & -{\int }_{{\Gamma }_{\sigma }}[{\hat{N}}_{nn}\delta {λ̶˜}_{u_{0n}}+{\hat{N}}_{ns}\delta {λ̶˜}_{u_{0s}}-{\hat{M}}_{nn}\frac{\mathrm{\partial }\left(\delta {λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }n}-{\hat{M}}_{ns}\frac{\mathrm{\partial }\left(\delta {λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }s}+{\hat{Q}}_n\left(\delta {λ̶˜}_{{\mathrm{w}}_0}\right)]\mathrm{d}s,\end{array}

where \left({\hat{N}}_{nn},{\hat{N}}_{ns}\right) and \left({\hat{M}}_{nn},{\hat{M}}_{ns}\right) are the in-plane forces resultants and moment resultants in the local coordinate \left(n,s,r\right). They are calculated by Eq. (15).

(15)\left\{\begin{array}{c}{\hat{N}}_{nn}\\ {\hat{N}}_{ns}\end{array}\right\}={\int }_{-h/2}^{h/2}\left\{\begin{array}{c}{\hat{\sigma }}_{nn}\\ {\hat{\sigma }}_{ns}\end{array}\right\}\mathrm{d} \mathit{\text{z}}, \left\{\begin{array}{c}{\hat{M}}_{nn}\\ {\hat{M}}_{ns}\end{array}\right\}={\int }_{-h/2}^{h/2}\left\{\begin{array}{c}{\hat{\sigma }}_{nn}\\ {\hat{\sigma }}_{ns}\end{array}\right\}\mathit{\text{z}}\mathrm{d}\mathit{\text{z}}.

And {\tau }_x and {\tau }_y are the direction if the unit normal vector \hat{n} is positioned at an angle {\tilde{\psi }}^{\ast } from x-axis and given in Eq. (16).

(16){\tau }_x=\mathrm{c}\mathrm{o}\mathrm{s}\left({\tilde{\psi }}^{\ast }\right), {\tau }_y=\mathrm{s}\mathrm{i}\mathrm{n}\left({\tilde{\psi }}^{\ast }\right).

The comma defined by subscripts in Eq. (14) represents differentiation concerning the subscripts {\stackrel{x}{N}}_{xx}=\mathrm{\partial }{N}_{xx}/\mathrm{\partial }x, {\stackrel{xx}{M}}_{xx}={\mathrm{\partial }}^2{M}_{xx}/\mathrm{\partial }{x}^2, and so on.

And the two terms of \mathrm{\aleph }\left({λ̶˜}_{{\mathrm{w}}_0}\right) and \rho \left({λ̶˜}_{{\mathrm{w}}_0}\right) are shown in Eq. (17).

(17)\begin{array}{rl}\mathrm{\aleph }\left({λ̶˜}_{{\mathrm{w}}_0}\right)=& \frac{\mathrm{\partial }}{\mathrm{\partial }x}\left(N_{xx}\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }x}+N_{xy}\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }y}\right)+\frac{\mathrm{\partial }}{\mathrm{\partial }y}\left(N_{xy}\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }x}+N_{yy}\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }y}\right),\\ \rho \left({λ̶˜}_{{\mathrm{w}}_0}\right)=& \left(N_{xx}\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }x}+N_{xy}\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }y}\right){\tau }_x+\left(N_{xy}\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }x}+N_{yy}\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }y}\right){\tau }_y,\end{array}

where (N_{xx},N_{yy} ,N_{xy} ) are the in-plan forces resultants and (M_{xx},M_{yy}, M_{xy} ) are the moment resultants. And they are calculated from Eqs. (18)–(21):

(18)\left\{\begin{array}{c}{N}_{xx}\\ N_{yy}\\ N_{xy}\end{array}\right\}=\sum _{k=1}^N{\int }_{{\mathit{\text{z}}}_k}^{{\mathit{\text{z}}}_k+1}\left\{\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\sigma }_{xy}\end{array}\right\}\mathrm{d}\mathit{\text{z}}= \sum _{k=1}^N {\int }_{{\mathit{\text{z}}}_k}^{{\mathit{\text{z}}}_k+1}{\left[\begin{array}{ccc}{\overline{R}}_{11}& {\overline{R}}_{12}& {\overline{R}}_{16}\\ {\overline{R}}_{12}& {\overline{R}}_{22}& {\overline{R}}_{26}\\ {\overline{R}}_{16}& {\overline{R}}_{26}& {\overline{R}}_{66}\end{array}\right]}^{\left(k\right)}\left\{\begin{array}{c}{\epsilon }_{xx}^{\left(0\right)}+\mathit{\text{z}}{\epsilon }_{xx}^{\left(1\right)}\\ {\epsilon }_{yy}^{\left(0\right)}+\mathit{\text{z}}{\epsilon }_{yy}^{\left(1\right)}\\ {\Gamma }_{yy}^{\left(0\right)}+\mathit{\text{z}}{\Gamma }_{yy}^{\left(1\right)}\end{array}\right\}\mathrm{d}\mathit{\text{z}},

(19)\left\{\begin{array}{c}{N}_{xx}\\ N_{yy}\\ N_{xy}\end{array}\right\}=\left[\begin{array}{ccc}{\overline{G}}_{11}& {\overline{G}}_{12}& {\overline{G}}_{16}\\ {\overline{G}}_{12}& {\overline{G}}_{22}& {\overline{G}}_{26}\\ {\overline{G}}_{16}& {\overline{G}}_{26}& {\overline{G}}_{66}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_{xx}^{\left(0\right)}\\ {\epsilon }_{yy}^{\left(0\right)}\\ {\Gamma }_{yy}^{\left(0\right)}\end{array}\right\}+\left[\begin{array}{ccc}{\overline{L}}_{11}& {\overline{L}}_{12}& {\overline{L}}_{16}\\ {\overline{L}}_{12}& {\overline{L}}_{22}& {\overline{L}}_{26}\\ {\overline{L}}_{16}& {\overline{L}}_{26}& {\overline{L}}_{66}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_{xx}^{\left(1\right)}\\ {\epsilon }_{yy}^{\left(1\right)}\\ {\Gamma }_{yy}^{\left(1\right)}\end{array}\right\},

And the moment resultants are given in Eqs. (20) and (21).

(20)$$\begin{gathered} \left\{\begin{array}{c}{M}_{xx}\\ M_{yy}\\ M_{xy}\end{array}\right\}=\sum _{k=1}^N{\int }_{{\mathit{\text{z}}}_k}^{{\mathit{\text{z}}}_k+1}\left\{\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}{\sigma }_{xy}\end{array}\right\}\mathit{\text{z}}\mathrm{d}\mathit{\text{z}} \\ =\sum _{k=1}^N{\int }_{{\mathit{\text{z}}}_k}^{{\mathit{\text{z}}}_k+1} {\left[\begin{array}{ccc}{\overline{R}}_{11}& {\overline{R}}_{12}& {\overline{R}}_{16}\\ {\overline{R}}_{12}& {\overline{R}}_{22}& {\overline{R}}_{26}\\ {\overline{R}}_{16}& {\overline{R}}_{26}& {\overline{R}}_{66}\end{array}\right]}^{\left(k\right)}\left\{\begin{array}{c}{\epsilon }_{xx}^{\left(0\right)}+\mathit{\text{z}}{\epsilon }_{xx}^{\left(1\right)}\\ {\epsilon }_{yy}^{\left(0\right)}+\mathit{\text{z}}{\epsilon }_{yy}^{\left(1\right)}\\ {\Gamma }_{yy}^{\left(0\right)}+\mathit{\text{z}}{\Gamma }_{yy}^{\left(1\right)}\end{array}\right\}\mathit{\text{z}}\mathrm{d} \mathit{\text{z}}, \end{gathered}$$

(21)\left\{\begin{array}{c}{M}_{xx}\\ M_{yy}\\ M_{xy}\end{array}\right\}=\left[\begin{array}{ccc}{\overline{L}}_{11}& {\overline{L}}_{12}& {\overline{L}}_{16}\\ {\overline{L}}_{12}& {\overline{L}}_{22}& {\overline{L}}_{26}\\ {\overline{L}}_{16}& {\overline{L}}_{26}& {\overline{L}}_{66}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_{xx}^{\left(0\right)}\\ {\epsilon }_{yy}^{\left(0\right)}\\ {\Gamma }_{yy}^{\left(0\right)}\end{array}\right\}+\left[\begin{array}{ccc}{\overline{J}}_{11}& {\overline{J}}_{12}& {\overline{J}}_{16}\\ {\overline{J}}_{12}& {\overline{J}}_{22}& {\overline{J}}_{26}\\ {\overline{J}}_{16}& {\overline{J}}_{26}& {\overline{J}}_{66}\end{array}\right]\left\{\begin{array}{c}{\epsilon }_{xx}^{\left(1\right)}\\ {\epsilon }_{yy}^{\left(1\right)}\\ {\Gamma }_{yy}^{\left(1\right)}\end{array}\right\},

where \left({\overline{G}}_{ij},{\overline{H}}_{ij},{\overline{J}}_{ij}\right) are defined by Eq. (22).

(22)\begin{array}{rl}& {\overline{G}}_{ij}=\sum _{k=1}^N{\overline{R}}_{ij}^{\left(k\right)}\left({\mathit{\text{z}}}_{k+1}-{\mathit{\text{z}}}_k\right),{\overline{L}}_{ij}=\frac{1}{2}\sum _{k=1}^N{\overline{R}}_{ij}^{\left(k\right)}\left({\mathit{\text{z}}}_{k+1}^2-{\mathit{\text{z}}}_k^2\right),\\ & {\overline{J}}_{ij}=\frac{1}{3}{\sum }_{k=1}^N{\overline{R}}_{ij}^{\left(k\right)}\left({\mathit{\text{z}}}_{k+1}^3-{\mathit{\text{z}}}_k^3\right).\end{array}

Equation (14) includes two terms. The whole domain {\mathrm{\Omega }}_0 is the first term and the second is along the boundary {\mathrm{\Gamma }}_{\sigma }. Thus, to satisfy Eq. (14), both terms must be zero. The Euler−Lagrange equations of CLPT as shown in Eq. (14) can be derived by setting the coefficients of \left(\delta {λ̶˜}_{{\mathrm{u}}_0},\delta {λ̶˜}_{{\mathrm{v}}_0},\delta {λ̶˜}_{{\mathrm{w}}_0}\right) to zeros separately as given Eq. (23).

(23)\begin{array}{rl}& \delta {λ̶˜}_{{\mathrm{u}}_0}:\frac{\mathrm{\partial }{N}_{xx}}{\mathrm{\partial }x}+\frac{\mathrm{\partial }{N}_{xy}}{\mathrm{\partial }y}=0,\delta {λ̶˜}_{{\mathrm{v}}_0}:\frac{\mathrm{\partial }{N}_{xy}}{\mathrm{\partial }x}+\frac{\mathrm{\partial }{N}_{yy}}{\mathrm{\partial }y}=0,\\ & \delta {λ̶˜}_{{\mathrm{w}}_0}:\frac{{\mathrm{\partial }}^2{M}_{xx}}{\mathrm{\partial }{x}^2}+2\frac{{\mathrm{\partial }}^2{M}_{xy}}{\mathrm{\partial }y\mathrm{\partial }x}+\frac{{\mathrm{\partial }}^2{M}_{yy}}{\mathrm{\partial }{y}^2}+\mathrm{\aleph }\left({λ̶˜}_{{\mathrm{w}}_0}\right)+q=0.\end{array}

By substituting the values of Eqs. (19) and (21) into Eq. (23), we can obtain a system of differential equations It can be written in a compact form as shown in Eqs. (24)–(26).

(24)\begin{array}{rl}& {\overline{G}}_{11}\left(\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{u}}_0}}{\mathrm{\partial }{x}^2}+\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^2}\right)+{\overline{G}}_{12}\left(\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{v}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}+\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }y}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}\right)+{\overline{G}}_{16}(\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{u}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}+\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{v}}_0}}{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^2}\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }y}\\ +& \frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }y})-{\overline{L}}_{11}\left(\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^3}\right)-{\overline{L}}_{12}\left(\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }{y}^2}\right)-2{\overline{L}}_{16}\left(\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^2\mathrm{\partial }y}\right)+{\overline{G}}_{16}\left(\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{u}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}+\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}\right)\\ +& {\overline{G}}_{26}\left(\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{v}}_0}}{\mathrm{\partial }{y}^2}+\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }y}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{y}^2}\right)+{\overline{G}}_{66}\left(\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{u}}_0}}{\mathrm{\partial }{y}^2}+\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{v}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}+\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }y}+\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{y}^2}\right)\\ -& {\overline{L}}_{16}\left(\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^2\mathrm{\partial }y}\right)-{\overline{L}}_{26}\left(\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{y}^3}\right)-2{\overline{L}}_{66}\left(\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }{y}^2}\right)=0,\end{array}

(25)\begin{array}{rl}& {\overline{G}}_{16}\left(\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{u}}_0}}{\mathrm{\partial }{x}^2}+\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^2}\right)+{\overline{G}}_{26}\left(\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{v}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}+\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }y}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}\right)+{\overline{G}}_{66}(\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{u}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}+\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{v}}_0}}{\mathrm{\partial }{x}^2}+\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^2}\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }y}\\ +& \frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }y})-{\overline{L}}_{16}\left(\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^3}\right)-{\overline{L}}_{26}\left(\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }{y}^2}\right)-2{\overline{L}}_{66}\left(\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^2\mathrm{\partial }y}\right)+{\overline{G}}_{12}\left(\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{u}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}+\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}\right)\\ +& {\overline{G}}_{22}\left(\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{v}}_0}}{\mathrm{\partial }{y}^2}+\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }y}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{y}^2}\right)+{\overline{G}}_{26}\left(\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{u}}_0}}{\mathrm{\partial }{y}^2}+\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{v}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}+\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }y}+\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{y}^2}\right)\\ -& {\overline{L}}_{12}\left(\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^2\mathrm{\partial }y}\right)-{\overline{L}}_{22}\left(\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{y}^3}\right)-2{\overline{L}}_{26}\left(\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }{y}^2}\right)=0,\end{array}

(26)\begin{array}{rl}& {\overline{L}}_{11}\left(\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{u}}_0}}{\mathrm{\partial }{x}^3}+\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^2}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^2}+\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x}\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^3}\right)+{\overline{L}}_{12}\left(\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{v}}_0}}{\mathrm{\partial }{x}^2\mathrm{\partial }y}+\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}+\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }y}\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^2\mathrm{\partial }y}\right)\\ +& {\overline{L}}_{16}\left(\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{u}}_0}}{\mathrm{\partial }{x}^2\mathrm{\partial }y}+\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{v}}_0}}{\mathrm{\partial }{x}^3}+\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^3}\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }y}+2\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^2}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}+\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x}\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^2\mathrm{\partial }y}\right)\\ -& {\overline{J}}_{11}\left(\frac{{\mathrm{\partial }}^4{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^4}\right)-{\overline{J}}_{12}\left(\frac{{\mathrm{\partial }}^4{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^2\mathrm{\partial }{y}^2}\right)-2{\overline{J}}_{16}\left(\frac{{\mathrm{\partial }}^4{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^3\mathrm{\partial }y}\right)+2{\overline{L}}_{16}\left(\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{u}}_0}}{\mathrm{\partial }{x}^2\mathrm{\partial }y}+\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^2}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}+\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x}\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^2\mathrm{\partial }y}\right)\\ +& 2{\overline{L}}_{26}\left(\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{v}}_0}}{\mathrm{\partial }x\mathrm{\partial }{y}^2}+\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{y}^2}+\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }y}\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }{y}^2}\right)+2{\overline{L}}_{66}(\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{u}}_0}}{\mathrm{\partial }x\mathrm{\partial }{y}^2}+\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{v}}_0}}{\mathrm{\partial }{x}^2\mathrm{\partial }y}+\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^2\mathrm{\partial }y}\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }y}\\ +& \frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}+\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^2}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{y}^2}+\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x}\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }{y}^2})-2{\overline{J}}_{16}\left(\frac{{\mathrm{\partial }}^4{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^3\mathrm{\partial }y}\right)-2{\overline{J}}_{26}\left(\frac{{\mathrm{\partial }}^4{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }{y}^3}\right)-4{\overline{J}}_{66}\left(\frac{{\mathrm{\partial }}^4{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^2\mathrm{\partial }{y}^2}\right)\\ +& {\overline{L}}_{12}\left(\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{u}}_0}}{\mathrm{\partial }x\mathrm{\partial }{y}^2}+\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}+\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x}\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }{y}^2}\right)+{\overline{L}}_{22}\left(\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{v}}_0}}{\mathrm{\partial }{y}^3}+\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{y}^2}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{y}^2}+\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }y}\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{y}^3}\right)\\ +& {\overline{L}}_{26}\left(\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{u}}_0}}{\mathrm{\partial }{y}^3}+\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{v}}_0}}{\mathrm{\partial }x\mathrm{\partial }{y}^2}+\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }{y}^2}\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }y}+2\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{y}^2}+\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x}\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{y}^3}\right)\\ -& {\overline{J}}_{12}\left(\frac{{\mathrm{\partial }}^4{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^2\mathrm{\partial }{y}^2}\right)-{\overline{J}}_{22}\left(\frac{{\mathrm{\partial }}^4{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{y}^4}\right)-2{\overline{J}}_{26}\left(\frac{{\mathrm{\partial }}^4{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }{y}^3}\right)+\mathrm{\aleph }\left({λ̶˜}_{{\mathrm{w}}_0}\right)+q=0,\end{array}

where \mathrm{\aleph }\left({λ̶˜}_{{\mathrm{w}}_0}\right) is defined in Eq. (17). It can be observed that the system of Eqs. (21)–(23) involves a set of variables to be determined, including \left({λ̶˜}_{{\mathrm{u}}_0},{λ̶˜}_{{\mathrm{v}}_0},{λ̶˜}_{{\mathrm{w}}_0}\right), where \left({λ̶˜}_{{\mathrm{u}}_0},{λ̶˜}_{{\mathrm{v}}_0},{λ̶˜}_{{\mathrm{w}}_0}\right) are mentioned in Eqs. (2)–(4).

They contain second-order spatial derivatives of \left({λ̶˜}_{{\mathrm{u}}_0},{λ̶˜}_{{\mathrm{v}}_0}\right) and fourth-order spatial derivatives of \left({λ̶˜}_{{\mathrm{w}}_0}\right), so the CLPT registers eighth-order theory. To solve the problem, we need to establish at least eight boundary conditions for the problem. The relationship between the global \left(x,y,\mathit{\text{z}}\right) and local coordinate \left(n,s,r\right) given in Eq. (27).

(27)\left\{\begin{array}{c}x\\ y\\ \mathit{\text{z}}\end{array}\right\}=\left[\begin{array}{ccc}\mathrm{c}\mathrm{o}\mathrm{s}{\tilde{\Psi }}^{\ast }& -\mathrm{s}\mathrm{i}\mathrm{n}{\tilde{\Psi }}^{\ast }& 0\\ \mathrm{s}\mathrm{i}\mathrm{n}{\tilde{\Psi }}^{\ast }& \mathrm{c}\mathrm{o}\mathrm{s}{\tilde{\Psi }}^{\ast }& 0\\ 0& 0& 1\end{array}\right]\left\{\begin{array}{c}n\\ s\\ r\end{array}\right\}.

Therefore, the displacements \left({λ̶˜}_{{\mathrm{u}}_{0n}},{λ̶˜}_{{\mathrm{u}}_{0s}}\right) are calculated from \left({λ̶˜}_{{\mathrm{u}}_0},{λ̶˜}_{{\mathrm{v}}_0}\right) as shown in Eq. (28).

(28){λ̶˜}_{{\mathrm{u}}_0}={\tau }_x{λ̶˜}_{{\mathrm{u}}_{0n}}-{\tau }_y{λ̶˜}_{{\mathrm{u}}_{0s}}, {λ̶˜}_{v_0}={\tau }_y{λ̶˜}_{{\mathrm{u}}_{0n}}+{\tau }_x{λ̶˜}_{{\mathrm{u}}_{0s}}.

The normal and tangential derivatives \left({λ̶˜}_{{\mathrm{w}}_{0,n}},{λ̶˜}_{{\mathrm{w}}_{0,s}}\right) are calculated by \left({λ̶˜}_{{\mathrm{w}}_{0,x}},{λ̶˜}_{{\mathrm{w}}_{0,y}}\right) given in Eq. (29).

(29)\begin{array}{rl}& \frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x}={\tau }_x\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }n}-{\tau }_y\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }s},\\ & \frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }y}={\tau }_y\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }n}+{\tau }_x\frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }s}.\end{array}

According to the above relations, we can rewrite the boundary in terms of \left({λ̶˜}_{{\mathrm{u}}_{0n}},{λ̶˜}_{{\mathrm{u}}_{0s}}\right) and \left({λ̶˜}_{{\mathrm{w}}_{0,n}},{λ̶˜}_{{\mathrm{w}}_{0,s}}\right) given in Eq. (30).

(30)\begin{array}{rl}& \left(N_{xx}{\tau }_x+N_{xy}{\tau }_y\right)\delta {λ̶˜}_{{\mathrm{u}}_0}+\left(N_{xy}{\tau }_x+N_{yy}{\tau }_u\right)\delta {λ̶˜}_{{\mathrm{v}}_0}\\ =& \left(N_{xx}{\tau }_x+N_{xy}{\tau }_y\right)\left({\tau }_x\delta {λ̶˜}_{{\mathrm{u}}_{0n}}-{\tau }_y\delta {λ̶˜}_{{\mathrm{u}}_{0s}}\right)\\ & +\left(N_{xy}{\tau }_x+N_{yy}{\tau }_y\right)\left({\tau }_y\delta {λ̶˜}_{{\mathrm{u}}_{0n}}+{\tau }_x\delta {λ̶˜}_{{\mathrm{u}}_{0s}}\right).\end{array}

The stresses under the two coordinates \left(x,y,\mathit{\text{z}}\right) and \left(n,s,r\right) are shown by Eq. (31)

(31)\left\{\begin{array}{c}{\sigma }_{nn}\\ {\sigma }_{ns}\end{array}\right\}=\left[\begin{array}{ccc}{\tau }_x^2& {\tau }_y^2& 2{\tau }_x{\tau }_y\\ -{\tau }_x{\tau }_y& {\tau }_x{\tau }_y& {\tau }_x^2-{\tau }_y^2\end{array}\right]\left\{\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\sigma }_{xy}\end{array}\right\}.

Hence, we have the relation of \left(N_{nn},N_{ns}\right), \left(M_{nn},M_{ns}\right) and \left(N_{xx},N_{yy},N_{xy}\right), \left(M_{xx},M_{yy},M_{xy}\right) as show in Eq. (32).

(32)\begin{array}{rl}& \left\{\begin{array}{c}{N}_{nn}\\ N_{ns}\end{array}\right\}=\left[\begin{array}{ccc}{\tau }_x^2& {\tau }_y^2& 2{\tau }_x{\tau }_y\\ -{\tau }_x{\tau }_y& {\tau }_x{\tau }_y& {\tau }_x^2-{\tau }_y^2\end{array}\right]\left\{\begin{array}{c}{N}_{xx}\\ N_{yy}\\ N_{xy}\end{array}\right\},\\ & \left\{\begin{array}{c}{M}_{nn}\\ M_{ns}\end{array}\right\}=\left[\begin{array}{ccc}{\tau }_x^2& {\tau }_y^2& 2{\tau }_x{\tau }_y\\ -{\tau }_x{\tau }_y& {\tau }_x{\tau }_y& {\tau }_x^2-{\tau }_y^2\end{array}\right]\left\{\begin{array}{c}{M}_{xx}\\ M_{yy}\\ M_{xy}\end{array}\right\}.\end{array}

Based on the above equations, the boundary integrals in Eq. (14) can be written as Eq. (33)

(33)\begin{array}{rl}0=& {\int }_{{\Gamma }_{\sigma }}[\left(N_{nn}-{\hat{M}}_{nn}\right)\delta {λ̶˜}_{{\mathrm{u}}_{0n}}+\left(N_{ns}-{\hat{N}}_{ns}\right)\delta {λ̶˜}_{{\mathrm{u}}_{0s}}\\ & +\left({\stackrel{x}{M}}_{xx}{\tau }_x+{\stackrel{y}{M}}_{xy}{\tau }_x+{\stackrel{y}{M}}_{yy}{\tau }_y+{\stackrel{x}{M}}_{xy}{\tau }_y+\rho \left({λ̶˜}_{{\mathrm{w}}_0}\right)\right)\delta {λ̶˜}_{{\mathrm{w}}_0}\\ & -\left(M_{nn}-{\hat{M}}_{nn}\right)\frac{\mathrm{\partial }\left(\delta {λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }n}-\left(M_{ns}-{\hat{M}}_{ns}\right)\frac{\mathrm{\partial }\left(\delta {λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }s}]\mathrm{d}s.\end{array}

Based on the condition of Eq. (33), The boundary conditions are summarized as the primary variables (displacement) are \left({λ̶¯}_{{\mathrm{u}}_n},{λ̶¯}_{{\mathrm{u}}_s},{λ̶˜}_{{\mathrm{w}}_0},{\displaystyle \frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }n}},{\displaystyle \frac{\mathrm{\partial }{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }s}}\right) and the second variables (forces) are \left(N_{nn},N_{ns},Q_n,M_{nn},M_{ns}\right), where Q_n is calculated by Eq. (34).

(34)Q_n= \left({\stackrel{x}{M}}_{xx}+{\stackrel{y}{M}}_{xy}\right){\tau }_x+ \left({\stackrel{y}{M}}_{yy}+{\stackrel{y}{M}}_{xy}\right){\tau }_y+ \rho \left({λ̶˜}_{{\mathrm{w}}_0}\right).

Equation (23) has a total spatial differential order of eight (second-order of {λ̶˜}_{{\mathrm{u}}_0}, second-order of {λ̶˜}_{{\mathrm{v}}_0}, and four-order of {λ̶˜}_{{\mathrm{w}}_0}). This implies that there should be only eight equations of boundary, however, here we have ten equations (five for variables of displacement and five for variables of force). To fix the inconsistency, integration by parts is applied given in Eq. (35).

(35)-\propto {\Gamma }_{\sigma }\left(M_{ns}\frac{\mathrm{\partial }\left(\delta {λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }s}\right)=\propto {\Gamma }_{\sigma }\left(\frac{\mathrm{\partial }{M}_{ns}}{\mathrm{\partial }s}\right)\mathrm{\partial }\left(\delta {λ̶˜}_{{\mathrm{w}}_0}\right)\mathrm{d}s- {\left[M_{ns}\delta {λ̶˜}_{{\mathrm{w}}_0}\right]}_{{\Gamma }_{\sigma }}.

When the stress boundary is closed {\left[M_{ns}\delta {λ̶˜}_{{\mathrm{w}}_0}\right]}_{{\Gamma }_{\sigma }}=0. This term now must be added to Q_n as shown in Eq. (36).

(36)V_n\equiv Q_n+ \frac{\mathrm{\partial }{M}_{ns}}{\mathrm{\partial }s}.

Which should be balanced by applied force {\hat{Q}}_n. So, the boundary of CLPT is given in Eq. (37).

(37)\left({λ̶¯}_{{\mathrm{u}}_n},{λ̶¯}_{{\mathrm{u}}_s},{λ̶˜}_{{\mathrm{w}}_0},\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }n}\right),\left(N_{nn},N_{ns},V_n,M_{nn}\right).

3 Physics informed neural networks for analysis of static composite laminates

3.1 Technical background

First, to understand how deep learning networks can solve differential equations, let us consider solving a problem in the form of the following differential equation given in Eq. (38).

(38)\{\begin{array}{rlr}& D\left(x,\frac{\mathrm{\partial }\left(\mathrm{\wp }\right)}{\mathrm{\partial }x},\frac{{\mathrm{\partial }}^2\left(\mathrm{\wp }\right)}{\mathrm{\partial }{x}^2},...,v\right)=0,& \mathrm{\forall }x\in \mathrm{\Omega },\\ & \mathrm{\wp }\left(x_0\right)=g\left(x_0\right),& \mathrm{\forall }{x}_0\in \mathrm{\partial }{\mathrm{\Omega }}_{\mathrm{\wp }},\end{array}

where D is a differential operator, \mathrm{\wp }\left(x\right) is the solution that we need to find in order to satisfy the differential equation and \nu is a vector of initial parameters. Specifically, the value of the function \mathrm{\wp }\left(x\right) must satisfy a predetermined condition for the values of \mathrm{\forall } x_0\in \mathrm{\partial }{\mathrm{\Omega }}_{\mathrm{\wp }}, where \mathrm{\partial } {\mathrm{\Omega }}_{\mathrm{\wp }} is a subset of \mathrm{\Omega }. g\left(x\right) is a known function. The equation above can be solved using analytical methods. However, a recently developed method can also solve this problem by seeking an alternative model. If the substitute model uses a DNN, this method is called PINN.

3.2 Fully-connected neural networks (FC-NNs)

PINN are a type of model that can include the underlying principles that govern a given data set, particularly those expressed by PDEs [1]. In essence, PINN are flexible function approximators that utilize an understanding of physical laws to enhance the learning process. NNs are a collection of algorithms based on the structure of biological brain networks, created to perform tasks like classification and regression. These networks come in various types, each with unique neuron connection patterns and architectures, such as FC-NNs. An FC-NN defines a mapping from the input layer to the output. In mathematical terms, every neuron is triggered by an activation function, which is expressed as \sigma \left(\right). The activation function is applied to a linearly transformed combination of the outputs from neurons in the previous layer. This transformation is achieved using a weight matrix W and a bias vector b which adjusts the input values before passing them through the activation function. Equation (39) describes the process of calculating the output values based on the parameters of the NN.

(39)

where K is the number of layer NN. Here if we consider the input data as the point coordinates \left(x,y\right), then predicted values of displacements \left({λ̶˜}_{{\mathrm{u}}_0},{λ̶˜}_{{\mathrm{v}}_0},{λ̶˜}_{{\mathrm{w}}_0}\right) to be calculated as the outputs of model as shown in Fig.4.

To simplify the representation of the mapping from FC-NNs shown in Fig.4, this process can be replaced by a surrogate model f_{\theta }\left(x,y\right) with the subscript \theta =\left[\left\{W\right\},\left\{b\right\}\right] is a collection of matrices of hyper and \left(x,y\right) is the coordinate of any point of laminate. So, f_{\theta }\left(x,y\right) can be defined by mapping f_{\theta }\left(x,y\right):=R^2\to R^3.

The model for predicting the displacement of \left({λ̶˜}_{{\mathrm{u}}_0},{λ̶˜}_{{\mathrm{v}}_0},{λ̶˜}_{{\mathrm{w}}_0}\right) can subsequently be defined iteratively as follows

Input layers: {\mathit{\text{z}}}_0= \left\{\begin{array}{c}x\\ y\end{array}\right\}\in R^2.

Hidden layers: ; R^2\to R^3.

Output layer: {\mathit{\text{z}}}^{\left(K\right)}= \left\{\begin{array}{c}{λ̶˜}_{{\mathrm{u}}_0}\\ {λ̶˜}_{{\mathrm{v}}_0}\\ {λ̶˜}_{{\mathrm{w}}_0}\end{array}\right\}=\left(W^{\left(K\right)}{\mathit{\text{z}}}^{\left(K-1\right)}+b^{\left(K\right)}\right) = ; R^2\to R^3.

The ability of NNs to approximate complex functions enables them to represent the full solution of PDEs effectively. This makes them powerful tools for solving a wide range of physics-based problems. It is important to note that the surrogate model’s predicted output values are continuously computed at each step by taking partial derivatives concerning the input variables. These derivatives play a crucial role in formulating the loss function, which ensures the accuracy and stability of the approximation. To construct this loss function, mathematical formulations such as strong or weak form equations are commonly employed. These equations, particularly the differential equations, serve as fundamental descriptors of the physical phenomena governing the problem. By embedding these governing equations into the learning process, PINN method ensures that the solution adheres to the underlying physical laws, making it a powerful and efficient approach for solving PDEs in various scientific and engineering applications.

3.3 Loss functions and implementation details

From the physical equations Eqs. (24)–(26) and (37), we can establish the objective function of the PINN model. The physical significance of the PINN model lies in its ability to fully describe and seek solutions with acceptable errors when applied to solve a system of differential equations. To find the best solution \theta = \left[\left\{W\right\},\left\{b\right\}\right] of f_{\theta }\left(x,y\right):=R^2\to R^3. PINN will be trained based on the minimum of the loss function. Notably, we can employ automatic differentiation to obtain derivatives appearing in PDE functions defined by Eqs. (24)– (26) and PDE for the boundary conditions based on Eq. (37). Automatic differentiation will evaluate the partial derivative of the output\left({λ̶˜}_{u_0},{λ̶˜}_{v_0},{λ̶˜}_{w_0}\right) concerning the input coordinates \left(x,y\right) through the surrogate model f_{\theta }\left(x,y\right):=R^2\to R^3. To complete the training process for PINN, we need to establish two loss functions. The first loss function is defined by the discrepancy of the system of three equations Eqs. (24)–(26). It is denoted as Los{s}_{\mathrm{P}\mathrm{D}\mathrm{E}\mathrm{s}}, and is computed as in Eq. (40).

(40)\begin{array}{rl}Los{s}_{\mathrm{P}\mathrm{D}\mathrm{E}}=& \frac{1}{Q_{\mathrm{c}}}{\sum }_1^{Q_{\mathrm{c}}}[{\Vert \left(\mathrm{E}\mathrm{q}.(24)\right)-0\Vert }^2\\ & +{\Vert (\mathrm{E}\mathrm{q}.(25))-0\Vert }^2+{\Vert \left(\mathrm{E}\mathrm{q}.(26)\right)-0\Vert }^2],\end{array}

where Q_{\mathrm{c}} are the internal collocation points for the residuals in the domain {\mathrm{\Omega }}_{\mathrm{u}} used for training the DNN. Please note that these points will not appear on the boundary domain {\mathrm{\Omega }}_{\mathrm{b}\mathrm{c}}. These points will serve the purpose of increasing the constraint level for the DNN network to achieve the required accuracy level. The total number of these points needs to be sufficiently large and distributed throughout the entire interior domain of the plate {\mathrm{\Omega }}_{\mathrm{u}}. Because through these points, a predictive trend will be established that satisfies the differential equations Eqs. (24)–(26). The second loss function is defined by the discrepancy of the given boundary conditions equation Eq. (37). Please note that we have a total of 8 boundary conditions equations, including 4 equations for boundary conditions on displacement (BCd) and 4 equations for the boundary conditions on stress (BCs). These will be denoted as Los{s}_{\mathrm{B}\mathrm{C}\mathrm{d}} and Los{s}_{\mathrm{B}\mathrm{C}\mathrm{s}}. These values are computed using the Eqs. (41) and (42).

(41)Los{s}_{\mathrm{B}\mathrm{C}\mathrm{d}}=\frac{1}{Q_{\mathrm{p}}} \sum _{i=1}^{Q_{\mathrm{p}}}\left({\Vert \left({λ̶˜}_{{\mathrm{u}}_n}-{λ̶^}_{{\hat{\mathrm{u}}}_n}\right)\Vert }^2+{\Vert \left({λ̶˜}_{{\mathrm{u}}_s}-{λ̶^}_{{\hat{\mathrm{u}}}_s}\right)\Vert }^2+{\Vert \left({λ̶˜}_{{\mathrm{w}}_0}-{λ̶^}_{{\hat{\mathrm{w}}}_0}\right)\Vert }^2+{\Vert \left(\frac{\mathrm{\partial }\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }n}-\frac{\mathrm{\partial }\left({λ̶^}_{{\hat{\mathrm{w}}}_0}\right)}{\mathrm{\partial }n}\right)\Vert }^2\right).

(42)Los{s}_{\mathrm{B}\mathrm{C}\mathrm{s}}=\frac{1}{Q_{\mathrm{p}}} \sum _{i=1}^{Q_{\mathrm{p}}}\left({\Vert \left(N_{nn}-{\hat{N}}_{nn}\right)\Vert }^2+{\Vert \left(N_{ns}-{\hat{N}}_{ns}\right)\Vert }^2+{\Vert \left(V_n-{\hat{N}}_n\right)\Vert }^2+{\Vert \left(M_{nn}-{\hat{M}}_{nn}\right)\Vert }^2\right).

The hat “^” symbol indicates the applied thickness-integrated forces and moments. The symbols without superscripts represent the predicted values generated by the DNN. And Q_{\mathrm{p}} are boundary collocation points for the residuals in the domain of boundary. These points are used to train the DNN network to solve the boundary conditions equations given by Eqs. (41) and (42). Typically, these points are evenly distributed along the edges of the boundary with equal spacing. Fig.5 illustrates the internal collocation points and boundary collocation points used to solve PDE with different BCs, and BCd.

By combining Eqs. (40)–(42), the loss function of this problem as shown in Eq. (43).

(43)Loss=Los{s}_{\mathrm{P}\mathrm{D}\mathrm{E}}+\left({\lambda }_{\mathrm{B}\mathrm{C}\mathrm{d}}Los{s}_{\mathrm{B}\mathrm{C}\mathrm{d}}+{\lambda }_{\mathrm{B}\mathrm{C}\mathrm{s}}Los{s}_{\mathrm{B}\mathrm{C}\mathrm{s}}\right),

where \left({\lambda }_{\mathrm{B}\mathrm{C}\mathrm{d}},{\lambda }_{\mathrm{B}\mathrm{C}\mathrm{s}}\right) are the weights of the residuals on displacement boundary domain and stress boundary domain. Adding these two coefficients is intended to address the issue of unbalance in computing the derivatives of the loss functions Los{s}_{\mathrm{P}\mathrm{D}\mathrm{E}} and \left(Los{s}_{\mathrm{B}\mathrm{C}\mathrm{d}}+Los{s}_{\mathrm{B}\mathrm{C}\mathrm{s}}\right). The imbalance between the loss functions is a term used to describe the disparity or imbalance between the components of the loss function in a machine learning model. In a model, loss functions are typically used to measure the difference between predicted values and the actual values of the data. When the loss functions have significant differences in terms of weight or importance in the model, an imbalance can occur, causing the model to overly focus on one part of the data and neglect the rest. This can lead to the model not performing well for parts of the data that are not accurately evaluated in the loss function [51].

It is very common to see a PINN trained using a data set of collocation points on the inside domain and the boundary domain. Therefore, the choice of these types of data plays an important role in successfully conducting the training process. Typically, there are two strategies to generate this data set: 1) the collocation points are generated initially and kept unchanged throughout the training process; 2) the data points are randomly sampled at each training epoch, meaning the data will change during the training process. For each specific problem, either strategy 1 or strategy 2 is used. In this study, the authors employed the first strategy to create the training data set. To introduce randomness into the data, the Mini-batch size technique is utilized, with its quantity being randomly selected for each training epoch. The training process terminates when the following equation is satisfied.

(44){\left[\left\{W\right\},\left\{b\right\}\right]}_{\text{best}}=\mathrm{a}\mathrm{r}\mathrm{g}\underset{W,b\in \mathbb{R}}{\mathrm{m}\mathrm{i}\mathrm{n}} Loss\left(\left\{x,y\right\},f_{\theta }\left(x,y\right),H_{\theta }\right),

where Loss is the loss function defined by Eq. (43), H_{\theta } is the learning algorithm. In this way, the best solution of \theta = \left[\left\{W\right\},\left\{b\right\}\right] will be exploited. It can be observed that H_{\theta } this is one of the most important factors of PINN as it determines the learning behavior of PINN to find the best model. In this study, the Adaptive Moment Estimation (Adam) algorithm [52] is used to update the hyperparameters after each epoch. It employs a parameter update technique that shares similarities with root mean square propagation but incorporates an additional momentum term. This algorithm can be illustrated in Eq. (45).

(45)\begin{array}{rl}& m{p}_k={\beta }_1\left(m{p}_{k-1}\right)+\left(1-{\beta }_1\right)\mathrm{\nabla }Loss\left({\theta }_k\right),\\ & \nu p_k={\beta }_2\left(v{p}_{k-1}\right)+\left(1-{\beta }_2\right){\left[\mathrm{\nabla }Loss\left({\theta }_k\right)\right]}^2,\\ & {\theta }_{k+1}={\theta }_k-\frac{\alpha \left(m{p}_k\right)}{\sqrt{\left(\nu p_k\right)}+\epsilon },\end{array}

where {\beta }_1 and {\beta }_2 are the gradient decay and squared gradient decay factor, respectively. \alpha is the learning rate.

For the surrogate model f_{\theta }\left(x,y\right):=R^2\to R^3, to construct a good model, there are many factors to consider, such as the number of hidden layers (and the corresponding number of nodes in each hidden layer), the selection of activation function, the number of collocation points, and so on. Regarding PINN, the literature contains a multitude of studies examining the impact of these hyperparameters [53–55]. By referring to previous studies, in this research, the network structure of the surrogate model f_{\theta }\left(x,y\right):=R^2\to R^3 is determined in Tab.1. And the process of finding the best surrogate model is implemented in Fig.6.

4 Numerical examples

4.1 Example 1–especially orthotropic plates [0°/90°/0°]

To clarify the effectiveness of PINN in applying to laminate plate problems, this study provides two numerical examples ranging from simple to complex. The first example focuses on orthotropic plates, specifically considering material distribution [0°/90°/0°], which bending-stretching coupling coefficients {\overline{H}}_{ij} and bending twisting coefficients {\overline{J}}_{16}= {\overline{J}}_{26}= 0. The uniform distribution load on the plate surface has a constant value of q_{\mathrm{t}}= q and this problem is shown in Fig.7.

Based on the assumptions, the strong form equation of especially orthotropic laminates is given in Eq. (46).

(46){\overline{J}}_{11}\frac{{\mathrm{\partial }}^4\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }{x}^4}+2 \left({\overline{J}}_{12}+2{\overline{J}}_{66}\right)\frac{{\mathrm{\partial }}^4\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }{x}^2\mathrm{\partial }{y}^2}+{\overline{J}}_{22}\frac{{\mathrm{\partial }}^4\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }{y}^4}=q.

To demonstrate that PINN can solve this problem with an acceptable error, we will first find an analytical Solution using the Navier solution. The Navier method is based on expanding the displacement using a double trigonometric (Fourier) series with unknown parameters. For solving Eq. (46), we admit for especially orthotropic rectangular laminates, and the boundary conditions are satisfied by following form of the transverse deflection as shown in Eq. (47).

(47){λ̶¯}_{\mathrm{w}}\left(x,y\right)= \sum _{n=1}^{\mathrm{\infty }}\sum _{m=1}^{\mathrm{\infty }}{W}^{mn} \mathrm{s}\mathrm{i}\mathrm{n}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }y\right)\mathrm{s}\mathrm{i}\mathrm{n}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }x\right),

where {\widetilde{\mathrm{{\rm Y}}}}^{\ast }=m\text{π}/ a,\mathrm{a}\mathrm{n}\mathrm{d}{\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }=n\text{π} /b. W^{mn} are calculated using equations from Eqs. (48)– (51)

(48)q\left(x,y\right)= \sum _{n=1}^{\mathrm{\infty }}\sum _{m=1}^{\mathrm{\infty }}{Q}^{mn} \mathrm{s}\mathrm{i}\mathrm{n}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }y\right)\mathrm{s}\mathrm{i}\mathrm{n}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }x\right),

where Q^{mn} is defined by Eq. (49).

(49)Q^{mn} =\frac{4}{ab}{\int }_0^b{\int }_0^{a}q \left(x,y\right)\mathrm{s}\mathrm{i}\mathrm{n}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }y\right)\mathrm{s}\mathrm{i}\mathrm{n}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }x\right)\mathrm{d}x\mathrm{d}y.

Substitution of the expansions Eqs. (48) and (49) into Eq. (47) yields, we obtain the form as shown in Eq. (50).

(50)\begin{array}{rl}\sum _{n=1}^{\mathrm{\infty }}\sum _{m=1}^{\mathrm{\infty }}& \{-W^{mn}[{\overline{J}}_{11}{\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }\right)}^4+2\left({\overline{J}}_{12}+2{\overline{J}}_{66}\right){\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }\right)}^2{\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }\right)}^2\\ & +{\overline{J}}_{22}{\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }\right)}^4]+Q^{mn}\}\mathrm{s}\mathrm{i}\mathrm{n}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }y\right)\mathrm{s}\mathrm{i}\mathrm{n}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }x\right)=0,\end{array}

W^{mn} can be calculated by Eq. (51).

(51)W^{mn} =\frac{Q^{mn}}{d^{mn}},

where d^{mn} is given by Eq. (52).

(52)d^{mn} =\frac{{\text{π}}^4}{b^4}\left[{\overline{J}}_{11}{m}^4{s}^4+2\left({\overline{J}}_{12}+2{\overline{J}}_{66}\right)m^2{n}^2{s}^2+{\overline{J}}_{22}{n}^4\right],

where s=b/a. So, the solution of Eq. (46) becomes Eq. (53).

(53){λ̶˜}_{{\mathrm{w}}_0}=\sum _{n=1}^{\mathrm{\infty }}\sum _{m=1}^{\mathrm{\infty }}\frac{Q^{mn}}{d^{mn}}\mathrm{s}\mathrm{i}\mathrm{n}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }y\right)\mathrm{s}\mathrm{i}\mathrm{n}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }x\right).

The stresses at each kth lamina is calculated by Eq. (54)

(54)\begin{array}{rl}& {\left\{\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\sigma }_{xy}\end{array}\right\}}^{\left(k\right)}=-\mathit{\text{z}}\left[\begin{array}{ccc}{\overline{R}}_{11}& {\overline{R}}_{12}& 0\\ {\overline{R}}_{12}& {\overline{R}}_{22}& 0\\ 0& 0& {\overline{R}}_{66}\end{array}\right]\left\{\begin{array}{c}{\displaystyle \frac{{\mathrm{\partial }}^2\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }{x}^2}}\\ {\displaystyle \frac{{\mathrm{\partial }}^2\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }{y}^2}}\\ 2{\displaystyle \frac{{\mathrm{\partial }}^2\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }x\mathrm{\partial }y}}\end{array}\right\}=\mathit{\text{z}}\sum _{n=1}^{\mathrm{\infty }}\sum _{m=1}^{\mathrm{\infty }}{W}^{mn}\\ & \cdot \left\{\begin{array}{c}\left({{\overline{R}}_{11}}^{\left(k\right)}{\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }\right)}^2+{{\overline{R}}_{12}}^{\left(k\right)}{\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }\right)}^2\right)\mathrm{s}\mathrm{i}\mathrm{n}\left({\displaystyle \frac{n\text{π}y}{{\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }}}\right)\mathrm{s}\mathrm{i}\mathrm{n}\left({\displaystyle \frac{m\text{π}x}{{\widetilde{\mathrm{{\rm Y}}}}^{\ast }}}\right)\\ \left({{\overline{R}}_{12}}^{\left(k\right)}{\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }\right)}^2+{{\overline{R}}_{22}}^{\left(k\right)}{\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }\right)}^2\right)\mathrm{s}\mathrm{i}\mathrm{n}\left({\displaystyle \frac{n\text{π}y}{{\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }}}\right)\mathrm{s}\mathrm{i}\mathrm{n}\left({\displaystyle \frac{m\text{π}x}{{\widetilde{\mathrm{{\rm Y}}}}^{\ast }}}\right)\\ \left(-2{{\overline{R}}_{66}}^{\left(k\right)}{\widetilde{\mathrm{{\rm Y}}}}^{\ast }{\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }\right)\mathrm{c}\mathrm{o}\mathrm{s}\left({\displaystyle \frac{n\text{π}y}{{\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }}}\right)\mathrm{c}\mathrm{o}\mathrm{s}\left({\displaystyle \frac{m\text{π}x}{{\widetilde{\mathrm{{\rm Y}}}}^{\ast }}}\right)\end{array}\right\}.\end{array}

4.2 Comparison the results of the Navier solution and the physics informed neural networks

To assess the reliability and accuracy of the PINN solution, especially orthotropic plates [0°/90°/0°] will be examined as the first example. The surveyed plate has two equal-length sides a/b=1, the thickness of each lamina h_i= a/ 100=b/100. The lamina’s material characteristics are considered to be: E_1/E_2= 25; G_{12}=G_{13}=0.5E_2 and v_{12}= 0.25. Laminates under uniformly distributed load q=q_{\mathrm{t}}= E_2/1000. To address the imbalance between \left(Los{s}_{\mathrm{B}\mathrm{C}\mathrm{d}}+Los{s}_{\mathrm{B}\mathrm{C}\mathrm{s}}\right) and Los{s}_{\mathrm{P}\mathrm{D}\mathrm{E}} the coefficients for calculating the loss function based on the boundary conditions are adjusted {\lambda }_{\mathrm{B}\mathrm{C}\mathrm{d}}={\lambda }_{\mathrm{B}\mathrm{C}\mathrm{s}}=0.1. The PINN method will be employed to solve the physical equations Eq. (46) and PINN method will be employed to simultaneously satisfy the boundary conditions of the problem. The PINN is set up with initial parameters as shown in Tab.1. The process of searching for the best parameters of the PINN model will be conducted through iterations with a total number of epochs set to 15000. The convergence trend of the loss functions Los{s}_{\mathrm{P}\mathrm{D}\mathrm{E}}, Los{s}_{\mathrm{B}\mathrm{C}\mathrm{d}}Los{s}_{\mathrm{B}\mathrm{C}\mathrm{s}} and Loss of total is illustrated in Fig.8. After the investigation, the authors found that with the number of epochs set to 15000, the loss functions achieve the necessary convergence, and the PINN model provides results with an acceptable error. The final values of the loss function achieved are as follows Los{s}_{\mathrm{P}\mathrm{D}\mathrm{E}} = 1.0871 × 10⁻⁶, Los{s}_{\mathrm{B}\mathrm{C}\mathrm{d}}= 4.409 × 10⁻⁷ and Los{s}_{\mathrm{B}\mathrm{C}\mathrm{s}}= 1.1322 × 10⁻⁵. And Loss of total given in Eq. (43), Loss=1.285 × 10⁻⁵.

After successfully training the PINN model with a total of 15000 epochs, we will obtain the values of the hyperparameters \theta = \left[\left\{W_{\mathrm{o}\mathrm{p}\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}\mathrm{d}}\right\},\left\{b_{\mathrm{o}\mathrm{p}\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}\mathrm{d}}\right\}\right]. Based on these hyperparameters, we will create a surrogate model that will be used to solve the problem in this example. The analysis of this problem using PINN will be conducted in two steps. 1) Computing the displacement at any point on the plate (including boundary points) based on Eq. (55).

(55)

2) Then, based on the known displacements, utilizing the automatic differentiation method, we can calculate the higher-order partial derivatives of the displacement concerning the positional coordinates. Stress at kth of the laminate plate will be calculated based on Eq. (56).

(56){\left\{\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\sigma }_{xy}\end{array}\right\}}^{\left(k\right)}=- \mathit{\text{z}}\left[\begin{array}{ccc}{\overline{R}}_{11}& {\overline{R}}_{12}& 0\\ {\overline{R}}_{12}& {\overline{R}}_{22}& 0\\ 0& 0& {\overline{R}}_{66}\end{array}\right]\left\{\begin{array}{c}{\displaystyle \frac{{\mathrm{\partial }}^2\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }{x}^2}}\\ {\displaystyle \frac{{\mathrm{\partial }}^2\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }{y}^2}}\\ 2{\displaystyle \frac{{\mathrm{\partial }}^2\left({λ̶˜}_{{\mathrm{w}}_0}\right)}{\mathrm{\partial }x\mathrm{\partial }y}}\end{array}\right\}.

The comparison results between using the PINN model and the Navier solution method will be conducted through the maximum values of displacement and maximum values of stress, which are simplified to nondimensionalized. The maximum deflection at the middle plane \left(x,y,\mathit{\text{z}}\right)=\left({\displaystyle \frac{a}{2}},{\displaystyle \frac{b}{2}},0\right), and the maximum normal stresses that occur at \left(x,y,\mathit{\text{z}}\right)=\left({\displaystyle \frac{a}{2}},{\displaystyle \frac{b}{2}},{\displaystyle \frac{h}{2}}\right) are the selected quantities for comparison. The maximum comparison values will be converted to nondimensionalized as given in Eq. (57) and shown in Tab.2.

(57)\begin{array}{rl}& \widehat{w_{xy}}={λ̶¯}_{{\mathrm{w}}_{xy}}\left(\frac{a}{2},\frac{b}{2},0\right)\left(\frac{E_2{h}^3}{a^4{q}_0}\right)\times 10^2,\\ & {\overline{\sigma }}_{xx}={\sigma }_{xx}\left(\frac{a}{2},\frac{b}{2},\frac{h}{2}\right)\left(\frac{h^2}{a^2{q}_0}\right),\\ & {\overline{\sigma }}_{yy}={\sigma }_{yy}\left(\frac{a}{2},\frac{b}{2},\frac{h}{2}\right)\left(\frac{h^2}{a^2{q}_0}\right).\end{array}

Based on the results in Tab.2, it can be observed that the values obtained through PINN approximate the numerical solution of the Navier solution. The maximum error value between the two methods recorded is 0.0049 at point \left({\displaystyle \frac{a}{2}},{\displaystyle \frac{b}{2}},0\right). To provide a more comprehensive comparison between the two methods, the absolute nondimensionalized displacement value at the center line of the plate (0\le x\le a,y=b) is used, the shape of the displacement results is shown in Fig.9(a). It can be observed that PINN yields excellent results and closely matches the values obtained by the Navier solution. Fig.9(b) and Fig.9(c) illustrate the absolute displacement of the entire plate, including the displacement at the positions of all edges. It can be observed that both the PINN and the Navier solution methods yield very similar results. In particular, the displacement at the edge positions approximates zero, which is consistent with the boundary conditions of the problem. The error spectrum between the two methods is shown in Fig.9(d), where very small error values approximating 6 × 10⁻³ are recorded.

As shown in Fig.9, the displacement of the plate obtained by the PINN method closely matches the numerical solution, both within the plate domain and along the plate edges. However, to comprehensively evaluate the plate’s behavior, in addition to the displacement values, the corresponding stresses {\overline{\sigma }}_{xx} and {\overline{\sigma }}_{yy} at different layer depths of the plate are also extracted for comparison. These values aim to assess whether the PINN model satisfies the force conditions through the automatic differentiation process based on the displacement. Therefore, by applying this method, we can compute the strain at any point given coordinates \left(x,y\right). Based on this strain, the stress values along the thickness direction of the plate {\overline{\sigma }}_{xx}, and {\overline{\sigma }}_{yy} will be calculated at the coordinate points \left({\displaystyle \frac{a}{2}},{\displaystyle \frac{b}{2}},-{\displaystyle \frac{h}{2}}\le \mathit{\text{z}}\le {\displaystyle \frac{h}{2}}\right). The comparison values are shown in Fig.10. Accordingly, the results indicate that the stresses {\overline{\sigma }}_{xx} and {\overline{\sigma }}_{yy} are very close to the Navier solution.

This example has demonstrated the effectiveness and reliability of the PINN method for analyzing a simple laminated plate problem. From the achieved results, it can be observed that the trained PINN model is capable of predicting displacement and stress values with an acceptable level of error. However, this example is presented in a simplified form. To further validate the effectiveness of the proposed method, a more general example of a cross-ply laminate plate will be investigated in the next section.

4.3 Example 2—Navier solution for antisymmetric cross-ply laminated plates \left[0^{\circ }/90^{\circ }/0^{\circ }/90^{\circ }\right]

In this example, a general asymmetric plate problem is studied with the material distribution \left[0^{\circ }/90^{\circ }/0^{\circ }/90^{\circ }\right]. This example is more complex and demonstrates a higher level of generality compared to Example 1. The investigated plate has equal lengths a/b=1, the thickness of each lamina in the plate. The lamina’s material characteristics are considered to be: E_1/E_2= 25; G_{12}=G_{13}=0.5E_2 and v_{12}= 0.25. Laminates under uniformly distributed load q=q_{\mathrm{t}}= E_2/1000. The geometry and boundary conditions of this example are shown in Fig.11.

The linear equations of motion of this example can be obtained from Eqs. (24)–(26) by setting the nonlinear terms to zero. We will obtain three reduced equations Eqs. (58)–(60).

(58)\begin{array}{rl}& {\overline{G}}_{11}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{u}}_0}}{\mathrm{\partial }{x}^2}+2{\overline{G}}_{16}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{u}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}+{\overline{G}}_{66}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{u}}_0}}{\mathrm{\partial }{y}^2}+{\overline{G}}_{16}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{v}}_0}}{\mathrm{\partial }{x}^2}+\left({\overline{G}}_{12}+{\overline{G}}_{66}\right)\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{v}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}+{\overline{G}}_{26}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{v}}_0}}{\mathrm{\partial }{y}^2}\\ -& {\overline{L}}_{11}\left[\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^3}+3{\overline{L}}_{16}\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^2\mathrm{\partial }y}+\left({\overline{L}}_{12}+2{\overline{L}}_{66}\right)\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }{y}^2}+{\overline{L}}_{26}\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{y}^3}\right]=0.\end{array}

(59)\begin{array}{rl}& {\overline{G}}_{16}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{u}}_0}}{\mathrm{\partial }{x}^2}+\left({\overline{G}}_{12}+{\overline{G}}_{66}\right)\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{u}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}+{\overline{G}}_{26}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{u}}_0}}{\mathrm{\partial }{y}^2}+{\overline{G}}_{66}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{v}}_0}}{\mathrm{\partial }{x}^2}+2{\overline{G}}_{66}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{v}}_0}}{\mathrm{\partial }x\mathrm{\partial }y}+{\overline{G}}_{22}\frac{{\mathrm{\partial }}^2{λ̶˜}_{{\mathrm{v}}_0}}{\mathrm{\partial }{y}^2}\\ -& {\overline{L}}_{16}\left[\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^3}+\left({\overline{L}}_{12}+2{\overline{B}}_{66}\right)\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^2\mathrm{\partial }y}+3{\overline{L}}_{26}\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }{y}^2}+{\overline{L}}_{22}\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{y}^3}\right]=0.\end{array}

(60)\begin{array}{rl}& {\overline{L}}_{11}\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{u}}_0}}{\mathrm{\partial }{x}^3}+3{\overline{L}}_{16}\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{u}}_0}}{\mathrm{\partial }{x}^2\mathrm{\partial }y}+\left({\overline{L}}_{12}+2{\overline{L}}_{66}\right)\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{u}}_0}}{\mathrm{\partial }x\mathrm{\partial }{y}^2}+{\overline{L}}_{26}\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{u}}_0}}{\mathrm{\partial }{y}^3}+{\overline{L}}_{16}\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{v}}_0}}{\mathrm{\partial }{x}^3}+\left({\overline{L}}_{12}+2{\overline{B}}_{66}\right)\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{v}}_0}}{\mathrm{\partial }{x}^2\mathrm{\partial }y}+3{\overline{L}}_{26}\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{v}}_0}}{\mathrm{\partial }x\mathrm{\partial }{y}^2}\\ +& {\overline{L}}_{22}\frac{{\mathrm{\partial }}^3{λ̶˜}_{{\mathrm{v}}_0}}{\mathrm{\partial }{y}^3}-{\overline{J}}_{11}\left[\frac{{\mathrm{\partial }}^4{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^4}+4{\overline{J}}_{16}\frac{{\mathrm{\partial }}^4{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^3\mathrm{\partial }y}+2\left({\overline{J}}_{12}+2{\overline{J}}_{66}\right)\frac{{\mathrm{\partial }}^4{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{x}^2\mathrm{\partial }{y}^2}+4{\overline{J}}_{26}\frac{{\mathrm{\partial }}^4{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }x\mathrm{\partial }{y}^3}+{\overline{J}}_{22}\frac{{\mathrm{\partial }}^4{λ̶˜}_{{\mathrm{w}}_0}}{\mathrm{\partial }{y}^4}\right]=0.\end{array}

To address the coupled system of PDEs along with the prescribed boundary conditions, analogous to the preceding example, the Navier solution will also be employed in this case. To fulfill the boundary constraints of the problem, the solutions of Eqs. (58)–(60) must be expressed in the form of a series expansion, as presented in Eq. (61).

(61)\begin{array}{rl}& {λ̶˜}_{{\mathrm{u}}_0}=\sum _{n=1}^{\mathrm{\infty }}\sum _{m=1}^{\mathrm{\infty }}{U}^{mn}\mathrm{s}\mathrm{i}\mathrm{n}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }y\right)\mathrm{c}\mathrm{o}\mathrm{s}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }x\right)\\ & {λ̶˜}_{{\mathrm{v}}_0}=\sum _{n=1}^{\mathrm{\infty }}\sum _{m=1}^{\mathrm{\infty }}{V}^{mn}\mathrm{c}\mathrm{o}\mathrm{s}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }y\right)\mathrm{s}\mathrm{i}\mathrm{n}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }x\right)\\ & {λ̶˜}_{{\mathrm{w}}_0}=\sum _{n=1}^{\mathrm{\infty }}\sum _{m=1}^{\mathrm{\infty }}{W}^{mn}\mathrm{s}\mathrm{i}\mathrm{n}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }y\right)\mathrm{s}\mathrm{i}\mathrm{n}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }x\right),\end{array}

where {\widetilde{\mathrm{{\rm Y}}}}^{\ast }=m\text{π}/ a,{\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }=n\text{π} /b and a set of (U^{mn},V^{mn},W^{mn} ) are the parameters to be determined based on the boundary conditions and by solving the algebraic equations given in Eq. (62)

(62)\left[\begin{array}{ccc}{e}_{11}& e_{12}& e_{13}\\ e_{12}& e_{22}& e_{23}\\ e_{13}& e_{23}& e_{33}\end{array}\right]\left\{\begin{array}{c}{U}^{mn}\\ V^{mn}\\ W^{mn}\end{array}\right\}=\left\{\begin{array}{c}0\\ 0\\ Q^{mn}\end{array}\right\},

where

(63)Q^{mn} =\frac{4}{ab}{\int }_0^b{\int }_0^{a}q \left(x,y\right)\mathrm{s}\mathrm{i}\mathrm{n}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }y\right)\mathrm{s}\mathrm{i}\mathrm{n}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }x\right)\mathrm{d}x\mathrm{d}y,

and e_{mn} is the coefficients whose explicit form is given in Eq. (64)

(64)\begin{array}{rl}& e_{11}={\overline{G}}_{11}{\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }\right)}^2+{\overline{G}}_{66}{\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }\right)}^2,\\ & e_{12}=\left({\overline{G}}_{12}+{\overline{G}}_{66}\right)\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }\right)\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }\right),\\ & e_{13}=-{\overline{H}}_{11}{\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }\right)}^3-({\overline{H}}_{12}+2{\overline{H}}_{66})\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }\right){\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }\right)}^2,\\ & e_{22}={\overline{G}}_{66}{\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }\right)}^2+{\overline{G}}_{22}{\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }\right)}^2,\\ & e_{23}=-\left({\overline{H}}_{12}+2{\overline{H}}_{66}\right){\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }\right)}^2\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }\right)-{\overline{H}}_{22}{\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }\right)}^3,\\ & e_{33}={\overline{J}}_{11}{\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }\right)}^4+2\left({\overline{J}}_{12}+2{\overline{J}}_{66}\right){\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }\right)}^2{\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }\right)}^2+{\overline{J}}_{22}{\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }\right)}^4.\end{array}

Based on Eq. (62), we can solve for the values of \left(U^{mn},V^{mn},W^{mn}\right). Then the final numerical solutions for displacements will be computed using Eq. (61). Once the displacement values have been computed, the stress values will be calculated using Eq. (65).

(65)\begin{array}{rl}{\left\{\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\sigma }_{xy}\end{array}\right\}}^{\left(k\right)}=& \sum _{m=1}^{\mathrm{\infty }}\sum _{n=1}^{\mathrm{\infty }}{\left[\begin{array}{ccc}{\overline{R}}_{11}& {\overline{R}}_{12}& 0\\ {\overline{R}}_{12}& {\overline{R}}_{22}& 0\\ 0& 0& {\overline{R}}_{62}\end{array}\right]}^{\left(k\right)}\\ & \left\{\begin{array}{c}\left(R_{mn}^{xx}+\mathit{\text{z}}{S}_{mn}^{xx}\right)\mathrm{s}\mathrm{i}\mathrm{n}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }y\right)\mathrm{s}\mathrm{i}\mathrm{n}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }x\right)\\ \left(R_{mn}^{yy}+\mathit{\text{z}}{S}_{mn}^{yy}\right)\mathrm{s}\mathrm{i}\mathrm{n}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }y\right)\mathrm{s}\mathrm{i}\mathrm{n}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }x\right)\\ \left(R_{mn}^{xy}+\mathit{\text{z}}{S}_{mn}^{xy}\right)\mathrm{c}\mathrm{o}\mathrm{s}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }y\right)\mathrm{c}\mathrm{o}\mathrm{s}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }x\right)\end{array}\right\},\end{array}

where \left(R_{mn}^{xx},R_{mn}^{yy},R_{mn}^{xy}\right) and \left(S_{mn}^{xx},S_{mn}^{yy},S_{mn}^{xy}\right) are calculated by Eq. (66).

(66)\begin{array}{rl}{\left\{\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\sigma }_{xy}\end{array}\right\}}^{\left(k\right)}=& \sum _{m=1}^{\mathrm{\infty }}\sum _{n=1}^{\mathrm{\infty }}{\left[\begin{array}{ccc}{\overline{R}}_{11}& {\overline{R}}_{12}& 0\\ {\overline{R}}_{12}& {\overline{R}}_{22}& 0\\ 0& 0& {\overline{R}}_{62}\end{array}\right]}^{\left(k\right)}\\ & \left\{\begin{array}{c}\left(R_{mn}^{xx}+\mathit{\text{z}}{S}_{mn}^{xx}\right)\mathrm{s}\mathrm{i}\mathrm{n}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }y\right)\mathrm{s}\mathrm{i}\mathrm{n}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }x\right)\\ \left(R_{mn}^{yy}+\mathit{\text{z}}{S}_{mn}^{yy}\right)\mathrm{s}\mathrm{i}\mathrm{n}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }y\right)\mathrm{s}\mathrm{i}\mathrm{n}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }x\right)\\ \left(R_{mn}^{xy}+\mathit{\text{z}}{S}_{mn}^{xy}\right)\mathrm{c}\mathrm{o}\mathrm{s}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast \ast }y\right)\mathrm{c}\mathrm{o}\mathrm{s}\left({\widetilde{\mathrm{{\rm Y}}}}^{\ast }x\right)\end{array}\right\}.\end{array}

4.4 Comparison the results of the Navier solution and the physics informed neural networks

The PINN method continues to be used for solving the problem in Example 2. The model parameters of PINN are set up as in Example 1. However, to address the unbalance between the loss functions \left(Los{s}_{\mathrm{B}\mathrm{C}\mathrm{d}}+Los{s}_{\mathrm{B}\mathrm{C}\mathrm{s}}\right) and Los{s}_{\mathrm{P}\mathrm{D}\mathrm{E}}, the authors adjust to achieve accuracy. It can be observed that due to the more complex boundary conditions in this example, the impact of the loss function associated with the BC \left(Los{s}_{\mathrm{B}\mathrm{C}\mathrm{d}}+Los{s}_{\mathrm{B}\mathrm{C}\mathrm{s}}\right) plays an important role in the optimization process of finding the best of \theta = \left[\left\{W_{\mathrm{o}\mathrm{p}\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}\mathrm{d}}\right\},\left\{b_{\mathrm{o}\mathrm{p}\mathrm{t}\mathrm{i}\mathrm{m}\mathrm{i}\mathrm{z}\mathrm{e}\mathrm{d}}\right\}\right]. According to the first example, with the selection of {\lambda }_{\mathrm{B}\mathrm{C}\mathrm{d}}={\lambda }_{\mathrm{B}\mathrm{C}\mathrm{s}}=0.1 and the results obtained are very promising. However, in this example, the complexity level of the problem was increased to achieve higher accuracy and precision of the model, accordingly, in this example, the coefficients \left({\lambda }_{\mathrm{B}\mathrm{C}\mathrm{d}},{\lambda }_{\mathrm{B}\mathrm{C}\mathrm{s}}\right) will be adjusted. Thus, the values of the two coefficients {\lambda }_{\mathrm{B}\mathrm{C}\mathrm{d}}={\lambda }_{\mathrm{B}\mathrm{C}\mathrm{s}}=0.1 if the number of epochs is less than 5000, and {\lambda }_{\mathrm{B}\mathrm{C}\mathrm{d}}={\lambda }_{\mathrm{B}\mathrm{C}\mathrm{s}}=1 if the number of epochs is greater than 5000. This is a novel proposal in this study based on the analysis process of the results with multiple independent runs. To compare the effectiveness of the proposed coefficients, the authors conducted the example with two PINN models, denoted as PINN (1) and PINN (2), with parameters set as specified in Tab.3, specifically for the two coefficients.

The purpose of this comparison was to assess the impact of different coefficient settings on the performance and outcomes of the PINN models. By varying these coefficients, the researchers aimed to evaluate their influence on the accuracy and reliability of the predictions made by the models.

Fig.12 illustrates the convergence trends of the loss functions as well as the influence of the coefficients {\lambda }_{\mathrm{B}\mathrm{C}\mathrm{d}} and {\lambda }_{\mathrm{B}\mathrm{C}\mathrm{s}}. It can be observed that for PINN (1) model at the number of Epoch = 5000, there is a clear division of the loss function values, the loss function of PINN (1), Los{s}_{\mathrm{B}\mathrm{C}\mathrm{s}} shows a significant improvement at epoch = 5000, where its increasing trend changes to a decreasing trend, which is sustained until the end of the epochs. Meanwhile, for the PINN (2) model, no significant improvement in the value of Los{s}_{\mathrm{B}\mathrm{C}\mathrm{s}} are observed throughout the entire training process. The result achieved by the PINN (2) model is Los{s}_{\mathrm{B}\mathrm{C}\mathrm{s}} = 7.6525 × 10⁻⁶ compared to the value of the PINN (1) model, which is Los{s}_{\mathrm{B}\mathrm{C}\mathrm{s}}= 2.3966 × 10⁻⁶. A similar trend is also observed in the value of Los{s}_{\mathrm{B}\mathrm{C}\mathrm{d}}. For the PINN (1) model at epoch = 5000, a sudden decrease is observed with a large magnitude. This favorable condition facilitates the achievement of a better value compared to the PINN (2) model. The final attained value of the PINN (1) model is Los{s}_{\mathrm{B}\mathrm{C}\mathrm{d}}=3.3033 × 10⁻⁶ comparing it to the value of the PINN (2) model, which is shown Los{s}_{\mathrm{B}\mathrm{C}\mathrm{d}}= 1.1861 × 10⁻⁴. Despite the PINN (1) model recording a value of Los{s}_{\mathrm{P}\mathrm{D}\mathrm{E}} = 1.7710 × 10⁻⁵ compared Los{s}_{\mathrm{P}\mathrm{D}\mathrm{E}} = 7.6248 × 10⁻⁶ in PINN (2) does not show any improvement. However, this value does not affect the total loss function calculated according to Eq. (43).

To clarify the effectiveness of adjusting the coefficients {\lambda }_{\mathrm{B}\mathrm{C}\mathrm{d}} and {\lambda }_{\mathrm{B}\mathrm{C}\mathrm{s}}, an addressing the imbalance between the loss functions dependent on the boundary conditions (Los{s}_{\mathrm{B}\mathrm{C}\mathrm{s}} and Los{s}_{\mathrm{B}\mathrm{C}\mathrm{s}}) compared to the loss function of the physical differential equation Los{s}_{\mathrm{P}\mathrm{D}\mathrm{E}}. L2 norm of the PINN (1) model and PINN (2) model are shown in Fig.13.

Similar to the observations in Fig.12, the L2 norm values of the boundary conditions functions (L2 norm BCs and L2 norm BCd) obtained from the PINN (1) model exhibit superior performance compared to those of the PINN (2) model. This suggests that PINN (1) more effectively enforces the given boundary constraints. In contrast, the PINN (2) model achieves a lower L2 norm value for the loss function associated with the differential equation, indicating its strength in satisfying the governing equation itself. However, it is important to note that if a model fails to adequately enforce the boundary conditions (BCs and BCd), the overall accuracy of the final solution may be compromised. The effectiveness of a PINN model, therefore, does not solely depend on minimizing the residuals of the PDE but also significantly on ensuring that the boundary conditions are properly satisfied. Neglecting the latter may result in a solution that deviates from practical expectations, ultimately reducing the reliability and applicability of the model. From this perspective, the PINN (1) model demonstrates better overall performance by maintaining a well-balanced approach to optimizing the loss functions of both the PDE and the boundary constraints. This balance is achieved through the incorporation of adjustable weighting coefficients {\lambda }_{\mathrm{B}\mathrm{C}\mathrm{d}} and {\lambda }_{\mathrm{B}\mathrm{C}\mathrm{s}} whose values, as shown in Tab.3, serve as a valuable reference for tuning PINN models. This approach is particularly beneficial for complex engineering applications, such as the analysis of laminate plates, where both governing equations and boundary conditions play crucial roles in achieving accurate and physically meaningful solutions.

The superior performance of the PINN (1) model over the PINN (2) model is clearly demonstrated through the nondimensionalized displacement results and the associated error levels when compared to the highly accurate Navier solution, as presented in Fig.14. A detailed examination reveals that the PINN (1) model is capable of generating displacement predictions that closely align with the Navier solution across the entire displacement field of the plate. This high degree of accuracy suggests that PINN (1) effectively captures the underlying physical behavior of the system. Conversely, the PINN (2) model exhibits significant discrepancies, failing to provide reliable displacement predictions for the plate. This inadequacy highlights the importance of properly addressing the imbalance between the two competing loss functions Los{s}_{\mathrm{B}\mathrm{C}\mathrm{d}}, Los{s}_{\mathrm{B}\mathrm{C}\mathrm{s}}, and Los{s}_{\mathrm{P}\mathrm{D}\mathrm{E}} in the model formulation. If these loss terms are not carefully calibrated, the optimization process may become biased, leading to inaccurate or even physically inconsistent solutions. Therefore, ensuring an appropriate balance between these loss functions is a critical aspect of successfully solving the problem, as it directly influences the overall accuracy and stability of the predicted displacement field. This insight underscores the necessity of developing systematic approaches for tuning loss weight parameters, particularly in complex engineering applications such as laminated plate analysis.

Similar to Example 1, the characteristic values to compare in this example are presented in Tab.4. Based on the results of the model with adjusted weights to address the imbalance in the loss function, the results of the PINN (1) model are used for comparison with the numerical solution. The comparative values include displacements and the maximum stress values. The results indicate that PINN (1) yields results very close to the Navier solution for the stress values {\overline{\sigma }}_{xx} and displacements {\overline{w}}_{xy}. The maximum absolute error is observed for the stress of {\overline{\sigma }}_{yy} with a magnitude of 0.1081.

The displacement distribution of the plate at the mid-plane, and the entire domain of the plate (including the boundary edges), and the error between the PINN method and the Navier solution are shown in Fig.15. The stress distribution at the point \left({\displaystyle \frac{a}{2}},{\displaystyle \frac{b}{2}},-{\displaystyle \frac{h}{2}}\le \mathit{\text{z}}\le {\displaystyle \frac{h}{2}}\right) is similar trend to Example 1, which is observed, at epochs = 15000. We can find a model that can predict the plate’s behavior with high accuracy, although in this example, the distribution and analysis are more complex compared to Example 1.

Fig.15 illustrates that the PINN (1) model yields results with an error margin that remains within an acceptable range when compared to the numerical solution for stress values. This observation substantiates the effectiveness of the proposed methodology, indicating that the approach developed in this study demonstrates a high degree of reliability in accurately capturing stress distributions.

5 Conclusions

This study introduces a novel computational framework utilizing the PINN approach for the static analysis of cross-ply composite laminates. By leveraging deep learning techniques in conjunction with the fundamental principles of solid mechanics, the proposed method provides an effective and reliable means of predicting stress and displacement distributions in composite laminates. The research approach involves formulating a system of PDEs based on static equilibrium principles and employing PINN to solve these equations in a data-efficient manner. To validate the effectiveness of the proposed methodology, two case studies with increasing complexity were analyzed. The results indicate that the PINN framework accurately captures the mechanical behavior of laminated plates, exhibiting high precision in predicting both stress and displacement fields. A key advantage of the approach is its ability to address the imbalance between different loss functions, ensuring stable and accurate convergence during the training process. The study highlights that optimizing the PINN model requires careful consideration of boundary conditions and an appropriate selection of weighting coefficients to balance the PDE residuals and boundary constraints effectively. The convergence trends of the loss functions further demonstrate the efficiency of the PINN method in reducing numerical errors, leading to solutions that closely match the reference Navier solutions. The observed maximum absolute error remains within a reasonable range, underscoring the robustness and reliability of the proposed approach for composite laminate analysis. Furthermore, visualizations of stress distributions and displacement patterns confirm the ability of the PINN framework to accurately capture the complex mechanical responses of laminated structures. Overall, this research contributes to the advancement of computational mechanics by introducing a physics-based NN framework tailored for composite materials. The proposed PINN methodology not only enhances the accuracy and efficiency of static analysis for laminated plates but also holds significant potential for addressing a broader range of engineering problems, including more complex structures and nonlinear material behaviors. Future research can further explore the adaptability of this approach to dynamic loading conditions, multi-scale modeling, and hybrid data-driven methodologies, expanding its applicability to real-world engineering challenges.

References

[1] Gibson R F. A review of recent research on mechanics of multifunctional composite materials and structures. Composite Structures, 2010, 92(12): 2793–2810

[2] Xiao Y, Qiao W, Fukuda H, Hatta H. The effect of embedded devices on structural integrity of composite laminates. Composite Structures, 2016, 153: 21–29

[3] Akbar M, Curiel-Sosa J. Piezoelectric energy harvester composite under dynamic bending with implementation to aircraft wingbox structure. Composite Structures, 2016, 153: 193–203

[4] Saeed M U, Li B B, Chen Z F. Mechanical effects of microchannels on fiber-reinforced composite structure. Composite Structures, 2016, 154: 129–141

[5] Cai S, Mao Z, Wang Z, Yin M, Karniadakis G E. Physics-informed neural networks (PINN) for fluid mechanics: A review. Acta Mechanica Sinica, 2021, 37(12): 1727–1738

[6] Yesilce Y. Effect of axial force on the free vibration of Reddy-Bickford multi-span beam carrying multiple spring-mass systems. Journal of Vibration and Control, 2010, 16(1): 11–32

[7] Karama M, Afaq K, Mistou S. Mechanical behaviour of laminated composite beam by the new multi-layered laminated composite structures model with transverse shear stress continuity. International Journal of Solids and Structures, 2003, 40(6): 1525–1546

[8] Saracoglu G, Ceritbinmez F, Özkan V, Yapici A. The mechanical effect of functionalization in laminated composite material. Aircraft Engineering and Aerospace Technology, 2023, 95(7): 1005–1015

[9] Garg N, Donough M J, Song C, Phillips A W, Prusty B G. Three-dimensional scaled boundary finite element method to study interfacial imperfections in thick laminated composite plates undergoing bi-axial bending. Engineering Analysis with Boundary Elements, 2023, 146: 34–43

[10] Dural E, Oyar F. Effect of delamination size, location and boundary conditions on the behavior of a laminated glass plate. Structures, 2023, 47: 121–133

[11] Minh H L, Sang-To T, Khatir S, Wahab M A, Gandomi A H, Cuong-Le T. Augmented deep neural network architecture for assessing damage severity in 3D concrete buildings under temperature fluctuations based on K-means optimization. Structures, 2023, 57: 105278

[12] Minh H L, Sang-To T, Abdel Wahab M, Cuong-Le T. A new metaheuristic optimization based on K-means clustering algorithm and its application to structural damage identification. Knowledge-Based Systems, 2022, 251: 109189

[13] Liu X, Qin J, Zhao K, Featherston C A, Kennedy D, Jing Y, Yang G. Design optimization of laminated composite structures using artificial neural network and genetic algorithm. Composite Structures, 2023, 305: 116500

[14] Hornik K, Stinchcombe M, White H. Multilayer feedforward networks are universal approximators. Neural Networks, 1989, 2(5): 359–366

[15] Chandrasekhar A, Suresh K. TOuNN: Topology optimization using neural networks. Structural and Multidisciplinary Optimization, 2021, 63(3): 1135–1149

[16] Minh H L, Sang-To T, Theraulaz G, Abdel Wahab M, Cuong-Le T. Termite life cycle optimizer. Expert Systems with Applications, 2023, 213: 119211

[17] Sang-To T, Le-Minh H, Abdel Wahab M, Thanh C L. A new metaheuristic algorithm: Shrimp and Goby association search algorithm and its application for damage identification in large-scale and complex structures. Advances in Engineering Software, 2023, 176: 103363

[18] Li Y, Minh H L, Cao M S, Qian X, Abdel Wahab M. An integrated surrogate model-driven and improved termite life cycle optimizer for damage identification in dams. Mechanical Systems and Signal Processing, 2024, 208: 110986

[19] Raissi M, Perdikaris P, Karniadakis G E. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational Physics, 2019, 378: 686–707

[20] ArendtWUrban K. Partial Differential Equations: An Introduction to Analytical and Numerical Methods. Berlin: Springer, 2023

[21] FlavinJ NRionero S. Qualitative Estimates for Partial Differential Equations: An Introduction. Boca Raton: CRC Press, 2020

[22] Li W, Bazant M Z, Zhu J. A physics-guided neural network framework for elastic plates: Comparison of governing equations-based and energy-based approaches. Computer Methods in Applied Mechanics and Engineering, 2021, 383: 113933

[23] Karniadakis G E, Kevrekidis I G, Lu L, Perdikaris P, Wang S, Yang L. Physics-informed machine learning. Nature Reviews. Physics, 2021, 3(6): 422–440

[24] Berg J, Nyström K. A unified deep artificial neural network approach to partial differential equations in complex geometries. Neurocomputing, 2018, 317: 28–41

[25] Wang Y, Sun J, Li W, Lu Z, Liu Y. CENN: Conservative energy method based on neural networks with subdomains for solving variational problems involving heterogeneous and complex geometries. Computer Methods in Applied Mechanics and Engineering, 2022, 400: 115491

[26] Jagtap A D, Kharazmi E, Karniadakis G E. Conservative physics-informed neural networks on discrete domains for conservation laws: Applications to forward and inverse problems. Computer Methods in Applied Mechanics and Engineering, 2020, 365: 113028

[27] Tang S, Feng X, Wu W, Xu H. Physics-informed neural networks combined with polynomial interpolation to solve nonlinear partial differential equations. Computers & Mathematics with Applications, 2023, 132: 48–62

[28] Gasick J, Qian X. Isogeometric neural networks: A new deep learning approach for solving parameterized partial differential equations. Computer Methods in Applied Mechanics and Engineering, 2023, 405: 115839

[29] Meng Z, Qian Q, Xu M, Yu B, Yıldız A R, Mirjalili S. PINN-FORM: A new physics-informed neural network for reliability analysis with partial differential equation. Computer Methods in Applied Mechanics and Engineering, 2023, 414: 116172

[30] Bai J, Rabczuk T, Gupta A, Alzubaidi L, Gu Y. A physics-informed neural network technique based on a modified loss function for computational 2D and 3D solid mechanics. Computational Mechanics, 2023, 71(3): 543–562

[31] De Florio M, Schiassi E, Furfaro R. Physics-informed neural networks and functional interpolation for stiff chemical kinetics. Chaos, 2022, 32(6): 063107

[32] Zhang P, Yin Z Y, Jin Y F, Yang J, Sheil B. Physics-informed multifidelity residual neural networks for hydromechanical modeling of granular soils and foundation considering internal erosion. Journal of Engineering Mechanics, 2022, 148(4): 04022015

[33] Cai S, Wang Z, Wang S, Perdikaris P, Karniadakis G E. Physics-informed neural networks for heat transfer problems. Journal of Heat Transfer, 2021, 143(6): 060801

[34] Bararnia H, Esmaeilpour M. On the application of physics informed neural networks (PINN) to solve boundary layer thermal-fluid problems. International Communications in Heat and Mass Transfer, 2022, 132: 105890

[35] Arzani A, Cassel K W, D’Souza R M. Theory-guided physics-informed neural networks for boundary layer problems with singular perturbation. Journal of Computational Physics, 2023, 473: 111768

[36] Almqvist A. Fundamentals of physics-informed neural networks applied to solve the Reynolds boundary value problem. Lubricants, 2021, 9(8): 82

[37] Sukumar N, Srivastava A. Exact imposition of boundary conditions with distance functions in physics-informed deep neural networks. Computer Methods in Applied Mechanics and Engineering, 2022, 389: 114333

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

[39] Zhuang X, Guo H, Alajlan N, Zhu H, Rabczuk T. Deep autoencoder based energy method for the bending, vibration, and buckling analysis of Kirchhoff plates with transfer learning. European Journal of Mechanics. A, Solids, 2021, 87: 104225

[40] Bastek J H, Kochmann D M. Physics-informed neural networks for shell structures. European Journal of Mechanics-A/Solids, 2023, 97: 104849

[41] Samaniego E, Anitescu C, Goswami S, Nguyen-Thanh V M, Guo H, Hamdia K, Zhuang X, Rabczuk T. An energy approach to the solution of partial differential equations in computational mechanics via machine learning: Concepts, implementation and applications. Computer Methods in Applied Mechanics and Engineering, 2020, 362: 112790

[42] Wang Y, Sun J, Rabczuk T, Liu Y. DCEM: A deep complementary energy method for linear elasticity. International Journal for Numerical Methods in Engineering, 2024, 125(24): e7585

[43] Peng L X, Sun J K, Tao Y P, Huang Z M. Bending analysis of thin plates with variable stiffness resting on elastic foundation via a two-network strategy physics-informed neural network method. Structures, 2024, 68: 107051

[44] Fallah A, Aghdam M M. Physics-informed neural network for bending and free vibration analysis of three-dimensional functionally graded porous beam resting on elastic foundation. Engineering with Computers, 2024, 40: 437–454

[45] Wang W, Thai H T. A physics-informed neural network framework for laminated composite plates under bending. Thin-Walled Structures, 2025, 210: 113014

[46] Vahab M, Haghighat E, Khaleghi M, Khalili N. A physics-informed neural network approach to solution and identification of biharmonic equations of elasticity. Journal of Engineering Mechanics, 2022, 148(2): 04021154

[47] Mishra B B, Kumar A, Samui P, Roshni T. Buckling of laminated composite skew plate using FEM and machine learning methods. Engineering Computations, 2021, 38(1): 501–528

[48] Furtado C, Pereira L F, Tavares R P, Salgado M, Otero F, Catalanotti G, Arteiro A, Bessa M A, Camanho P P. A methodology to generate design allowables of composite laminates using machine learning. International Journal of Solids and Structures, 2021, 233: 111095

[49] WangYSun JBaiJAnitescuCEshaghiM ZhuangXRabczuk TLiuY. Kolmogorov Arnold Informed neural network: A physics-informed deep learning framework for solving PDEs based on Kolmogorov Arnold Networks. 2024, arXiv: 2406.11045

[50] Reddy J N. A simple higher-order theory for laminated composite plates. Journal of Applied Mechanics, 1984, 51(4): 745–752

[51] Xiang Z, Peng W, Liu X, Yao W. Self-adaptive loss balanced Physics-informed neural networks. Neurocomputing, 2022, 496: 11–34

[52] KingmaD PBa J. Adam: A method for stochastic optimization. 2014, arXiv: 1412.6980

[53] Song J, Cao W, Liao F, Zhang W. VW-PINNs: A volume weighting method for PDE residuals in physics-informed neural networks. Acta Mechanica Sinica, 2024, 41: 324140

[54] Jeong H, Bai J, Batuwatta-Gamage C P, Rathnayaka C, Zhou Y, Gu Y T. A physics-informed neural network-based topology optimization (PINNTO) framework for structural optimization. Engineering Structures, 2023, 278: 115484

[55] Wu C, Zhu M, Tan Q, Kartha Y, Lu L. A comprehensive study of non-adaptive and residual-based adaptive sampling for physics-informed neural networks. Computer Methods in Applied Mechanics and Engineering, 2023, 403: 115671

RIGHTS & PERMISSIONS

Higher Education Press