A CNN-based surrogate model of isogeometric analysis in nonlocal flexoelectric problems

We proposed a convolutional neural network (CNN)-based surrogate model to predict the nonlocal response for flexoelectric structures with complex topologies. The input, i.e. the binary images, for the CNN is obtained by converting geometries into pixels, while the output comes from simulations of an isogeometric (IGA) flexoelectric model, which in turn exploits the higher-order continuity of the underlying non-uniform rational B-splines (NURBS) basis functions to fast computing of flexoelectric parameters, e.g., electric gradient, mechanical displacement, strain, and strain gradient. To generate the dataset of porous flexoelectric cantilevers, we developed a NURBS trimming technique based on the IGA model. As for CNN construction, the key factors were optimized based on the IGA dataset, including activation functions, dropout layers, and optimizers. Then the cross-validation was conducted to test the CNN’s generalization ability. Last but not least, the potential of the CNN performance has been explored under different model output sizes and the corresponding possible optimal model layout is proposed. The results can be instructive for studies on deep learning of other nonlocal mech-physical simulations.


Introduction
Unlike piezoelectricity coupling electric polarization directly with the strain, flexoelectricity is a nonlocal coupling of the polarization to the strain gradient, as shown in Fig. 1. While piezoelectricity exists only in non-centrosymmetric materials, flexoelectricity shows in nearly any material [27,53]. Promising applications of flexoelectricity have been emerging like energy harvesters absorbing the mechanical vibrations [22,32], actuators compatible for the semiconductor industries [8], self-powering microsensors [28,31], as well as domain tailoring and polarization switching [17,30]. In nano-scale systems, flexoelectricity is pronounced due to the arised large nonlocal strain gradients [2]. Simulation of flexoelectric materials is complicated due to the C 1 continuity requirement in the weak form. This numerical prerequisite complicates the implementation of the classical finite element methods based on Lagrange polynomials. Alternatively, the meshfree methods as presented in [3,40] or isogeometric analysis (IGA) [21] can be used which intrinsically fulfill the C 1 requirement.
However, simulating simple bulk materials is not yet enough to excavate the potentials of the nonlocal flexoelectricity. Possible mechanisms of enhanced flexoelectricity in dielectrics mainly includes inner strain, polar nano-regions, non-crystalline polar phases, residual ferroelectricity, and surface piezoelectricity [12,45]. Zhang et al. [55] found that the surface piezoelectricity can contribute about 70% of the enhanced flexoelectricity in BaTiO 3 ceramics. Moreover, the shape and configuration of materials have been reported affecting the piezoelectric and the flexoelectric effects [26,35,42]. To better develop flexoelectric devices based on the nonlocal electro-mechanical coupling effects, topology optimization considering the surface piezoelectricity has proved feasible for piezoelectric structures [34] and for flexoelectric structures [14,36]. As meshfree methods tend to be computationally more expensive, IGA have advantages in representing complex boundary representations [15] when adapting some improvements such as more sophisticated splines [37] or NURBS trimming technology [24].
Optimizing flexoelectricity with complex topologies required repetitive yet expensive simulations. Surrogate modeling serves right efficient for such problems, dramatically decreasing the computational cost by approximating the input-output relation of a system. Namely, given input parameters, e.g., initial/boundary/operational conditions, the quantities of interest (QoIs), such as electric potential, stress/ strain, and their integrals/gradients can be obtained rapidly without conducting the simulations. The existing surrogate modeling approaches can be roughly categorized into two classes: projection-based reduced order models (ROMs) and data-fit models [7]. Common surrogate modeling techniques include: (i) Response Surface Models [33], (ii) Kriging Methods [11,39], (iii) Radial Basis Functions [9,51], and (v) Support Vector Regression [6,46]. With the popularity of machine learning (ML), training surrogate models is also posed as a supervised learning problem [56], via Artificial Neural Networks (ANN) [5,13,52]. The Deep Neural Network (DNN), also called Deep learning, provides a better representation of non-linear input-output relationships compared to classical ANN [4]. Hamdia et al. [19] proved the accuracy of deep learning in solving flexoelectric problems. Goswami et al. [18] proposed a physics-informed DNN for modeling complex topologies like a fracture. Since first proposed in 1998 [29], the Convolutional Neural Network (CNN) has shown outstanding learning abilities in a wide range of applications including natural language [25,54], speech recognition [1,43], and computer vision [20,41].
Although ML models have also shown promise in many disciplines outside the realm of computer science (e.g., aquatic sciences [23], atmospheric science [38], and biomedical science [48]), the capability of CNN in extracting geometry-physics linkages, at present, is still difficult to satisfy the engineering applications, especially in the nonlocal flexoelectricity. Inspired by the exciting predicting performance of CNNs in other fields, in this paper, we developed a CNN-based surrogate model for the efficient analysis of flexoelectric structures with complex topology. Our method further demonstrated the underlying potential of optimizing flexoelectricity.
The paper is arranged as follows: Sect. 2 summarizes the boundary value problem (BVP) and associated IGA formulation. Section 3 describes the CNN model before we present several benchmark problems in Sect. 4. The manuscript finishes with some concluding remarks in Sect. 5.

Governing equation for flexoelectricity
The weak form of flexoelectricity for a linear dielectric solid can be derived from the electric enthalpy density h which is a function of the linear strain tensor ij , the electric field E i h ij , E i and the strain gradient jk,l [3]: where e ikl , ij and C ijkl denote the piezoelectricity tensor, the dielectric tensor and the elasticity tensor, respectively; ijkl = d ijkl − f ijkl is the flexoelectric tensor which includes the direct d ijkl and the reverse flexoelectricity tensor f ijkl , respectively. Integrating h over the domain , we obtain the total electric enthalpy where ̂ ij , ̃ ij and D i are defined as

IGA formulation
Our IGA is based on Non-Uniform rational B-splines (NURBS), which in turn originated from B-splines. The basis functions N i,p of p-order B-splines with (n + 1) control points are given as where is a knot vector in the parametric space. Introducing the weights i to B-splines base functions, the mechanical displacement and the electric potential on a p-order NURBS surface with (m + 1)(n + 1) control points can be defined as with The test functions have a similar structure. Substituting the trial functions and test functions into the weak form, Eq. (7) yields after some algebra to the final linear system of equations [15]: For details about IGA derivation, please refer to Appendix 1.

NURBS trimming technique
To facilitate the computing of flexoelectricity on complex topologies like inclusions and holes, we incorporated into our IGA formulation the NURBS trimming technique proposed by Kin et al. [24]. The flow of the trimming analysis is demonstrated in Fig. 2. First, the NURBS is meshed uu u = . in the parametric space, all the elements are identified and decomposed into triangular sub-elements. Then, for the numerical integration of the triangular sub-elements, the Jacobian matrix is calculated via the transformation of the sub-elements from the physical space to the Gaussian quadrature space. Finally, the element stiffness matrix and global stiffness matrix are assembled for solving the boundary value problem. Therefore, the key of the IGA on trimmed NURBS model lies in the identification of trimmed elements and the accuracy of the corresponding numerical integration.

Trimmed elements identification
In the NURBS geometries, the distance of a point to a trimming curve can be obtained by finding the smallest distance between the point to its projections on the curve, e.g., d 2 is minimum projection distance, and thus is the distance of the point P to the curve C(u), as shown in Fig. 3a.
Hence, the trimmed elements can be generally recognized by the distance relationship between the element's center, the element's vertices, and the curve, as shown in Fig. 3b. Possible conditions are as follows: (i) If d c < r in , the element is trimmed by the curve; (ii) If d c > r out , the element is normal; (iii) If r in ≤ d c ≤ r out , the element needs refinement with the Quadtree algorithm [24].

Trimmed elements decomposition
In Jacobian computing, the trimmed elements can be divided into 3 types, based on the number of vertices outside of the physical domain, as shown in Fig. 4. NEFEM (NURBS integral finite element method) [44] is used to process straight-sided triangular sub-elements and curved-sided triangular sub-elements, as shown in Appendix 2.

Case I: an infinite plate with a circular hole
Our IGA method with NURBS trimming was implemented in MATLAB. To validate the proposed approach, we first investigated a infinite plate with a circular hole of radius a = 1 m (case I), which is subjected to a subjected to a far field traction ( = 1 Pa ) in the x direction, as shown in Fig. 5. A finite portion of the infinite plate is considered for analysis and, due to the symmetry of the problem in Fig. 5a, only a quarter of the portion requires modelling in Fig. 5b. The plane stress conditions are assumed as our previous study [10], and, the elastic material properties used   Fig. 5b, the stresses and displacements for this benchmark are given in an analytical solution in [49] as and where G is the shear modulus and the Kolosov constant = (3 − )∕(1 − ) for the plane strain assumption. Boundary conditions specified in [49] are applied on the right and upper edges as shown in Fig. 5b.
The performance of the proposed IGA formulation is studied using meshes with 36, 100, 324 and 1156 nodes as shown in Fig. 6b-d and the NURBS basis is used for all meshes. The proposed method shows good convergence characteristics in this problem involving stress concentration. Figure 6e shows xx predicted along the left edge x = 0 showing comparable accuracy and smoothness, and close agreement with the analytical solution by Eq. (12), where xx ( = ∕2, r = a) = 3.

Case II: 2-D flexoelectric cantilever with multi-holes
Furthermore, the 2-D flexoelectric cantilever (case II) in Fig. 7 was investigated. The length L and the height h of the cantilever are 0.454 mm and 0.227 mm, respectively, As for case II, our MATLAB program solved the displacement and the electric potential in Eq. (11) using the material properties in Table 1. Both scenarios whether the beam in Fig. 7 is with holes or not, are simulated. The impact of holes on electric potential is shown in Fig. 8. More detailed interactions of the defects with the structure response were revealed by the contours of displacement, the strain, and the strain gradient in Figs. 9 and 10.

Dataset generation
Our CNN-based surrogate model is designed to predict the identical output from our IGA model, i.e., the basic unknowns including the displacement and the electric potential. The other physical fields like the strain/stress and the strain/stress gradient can be derived from the predicted basic unknowns. To simplify the generation of representative porous flexoelectric structures, the NURBS representation for circles is used to indicate the holes as shown in Fig. 11.
Similar to the boundary settings of Fig. 7, a cantilever is fixed at one end and applied with a point load of 100 mN at the top vertex of the other end, while the bottom edge is grounded. Then the MATLAB's built-in uniform distribution algorithm was utilized to generate the circles with arbitrary radius r in a defined range {r | 0.05h ≤ r ≤ 0.2h} and to guarantee an equal chance inside the defined limit boundary in Fig. 11. The corresponding material properties used are shown in Table 1. The physical meanings of these properties can be referred to in Eq. (1).
The dataset structure is shown in Fig. 12. For each sample, the CNN-based model input is a binary image recording the geometry information and the model output is an array of the nodal basic unknowns, i.e., the mechanical displacement components ( u 1 , u 2 ) and the electric potential ( ).
For the following study, we considered the number of holes (nb) ranging from 1 to 5. For each case, 200 samples were produced randomly. And three cases of model output numbers (with different h-refinements) were considered. In total, the dataset contains 3000 samples, as shown in Fig. 12. For most cases, the simulation is done in 2 min. But for samples with the high h-refinement, the calculation could take over 5 min due to the NURBS trimming analysis.

Data preprocessing
The dataset is randomly split into the training set (80% data) and the testing dataset (20% data) in the model training. For computational efficiency, the pixels of the binary images are reduced to 64 × 32 . The Min-Max normalization is applied to the physical fields data, which on the dataset E = {e i } is written as After the Min-Max normalization, a good balance between the model bias and the data variance is achieved, judging from the loss difference between test and train datasets in   Fig. 13. The metric for model loss adopted here is the mean square error (MSE), which is the sum of squared distances between targets e i and predictions ê i .

Model building
The proposed CNN-based surrogate model is shown in Fig. 14a. Realized in Python with PyTorch library, the model input is binary images of flexoelectric porous beams, while the output is the same as that from IGA simulations. When the training data are input in our model, three 2-D convolutional layers (Con2D) and three Max-Pooling2D layers are arranged alternately, to gradually extract the underline data features at different abstract extend as the network grows deeper, as shown in Fig. 14b. Then to avoid the common overfitting issue in CNN, the dropout layer is followed. In the end, the fully connected Dense layer is added according to make the learned features match the dimension of the output. Unlike classic pixel-to-pixel predictions in CNNs, our CNN model predicts nodal fields as the output array instead, which ameliorates the cost of reducing image resolution, so that the training efficiency and prediction accuracy is improved.

Activation function
ReLU (Rectified Linear Unit) activation functions exhibit certain biological meaning [16] and are widely used in deep neural network applications including image  recognition. The ReLU is used for all the layers in our models, except that the final hidden layer must use linear activation for prediction's purpose, as shown in Fig. 14c.

Optimizer
Under the same training settings (200 samples with 80% as training dataset, epochs = 50, batch size = 20), the Adaptive Moment Estimation (Adam) outperforms other common adaptive techniques, as shown in Table 2.

Dropout
Dropout is a regularization technique for neural network models to avoid overfitting [47]. Through dropout layers, random neurons are ignored during one training epoch, so their contribution to the downstream neurons and weight updates on the backward pass is temporally removed. Dropout layers are normally used after all convolutional layers and pooling layers as indicated in our CNN-based surrogate model (Fig. 15a). An alternative way is to drop out entire feature maps from the convolutional layer which are then not used during pooling. This is called spatial dropout [50], as illustrated in Fig. 15b.
A low dropout has a minimal effect while a high value may result in under-learning. In our model, the dropout value is set to 50%, which seems close to the optimal for a wide range of networks [47]. Under the same training settings (200 samples with 80% as training dataset, epochs = 50, batch size = 20), the validation loss of (a) is 0.0068, which is slightly better than that of (b), which is 0.0088. So our model adopted the normal dropout layers.

Results and analysis
Our CNN model could predict the identical physical unknowns as our IGA model at a stable accuracy by the model selection and hypermeter tuning. Details of training settings and results will be discussed in the following. We first tested the accuracy of our model by cross-validation, and then analyzed the potential of our model performance and possible optimal model layout.

Cross-validation
That k-fold cross-validation (CV) is a resampling procedure to estimate the generalization ability of the model on new data. Here k refers to the number of groups that the samples are to be split into. To validate our model, 10-fold crossvalidation was conducted, as shown in Fig. 16. Under the same training settings (200 samples with 80% as training

Model performance analysis
In building our CNN-based surrogate model, the prediction errors exist due to the discretization from the analytic NURBS into pixels, due to the reduction of the pixels in data preprocessing, and finally the CNN approximation. Meanwhile, with the specified CNN layers, our model may not be right the best for the given problem.
It would therefore be interesting to see how to minimize the total errors by adjusting the size of the model output and convolutional layers. The layers of our current model are summarized in Table 3. The size of the model output is determined by the aforementioned h-refinement. For a scalar physical field like electric potential, the output size ranges in {100, 324, 1156} while h is in {3, 4, 5} . The influence of   Fig. 17.
In the analysis, 1000 samples were randomly selected, where 80% samples were set as the training dataset, while the rest 20% as the testing dataset. For all cases, the MSE turned to converge after the first 100 epochs. The training time is comparable, around 200 seconds each. When h is 3, the model output size is 100, and the corresponding MSE is the smallest, 6 × 10 −4 . As the model output size increased, MSE jumped around 40 times bigger, which means the current model might have no potential for predicting under the larger output size. In fact, when h is 3 with 100 control points, the mesh is very coarse. Thus, the error from NURBS trimming of our IGA model is considerable compared to the finer cases.
To improve the performance of our basic CNN model in Table 3, the size of convolutional layers is our following concern. We denote the size of convolutional layers as scale = 1. There are about 4.4 × 10 5 hyperparameters for the basic model, and most are in the convolutional layers. In the following analysis, the neurons in all the three Conv2D layers are changed by the same scale. Considering the practical training time, the explored network scale is limited in {2 −4 , 2 −2 , 2 0 , 2 2 , 2 4 } . For a given scale 2 i , the number of hyperparameters in the corresponding model is about 4.4 ⋅ 2 i × 10 5 . 200 samples were randomly selected, where 80% samples were set as the training dataset, while the rest 20% as the testing dataset.
The MSE is evaluated for the first 100 epochs. For varying output sizes, the potential optimal size of convolutional layers differs, as shown in Table 4. While the original model scale is most suitable for the case that output size is 100, the increased model scale did make our model predict better in larger output sizes of 324 and 1156. The number of hyperparameters reflects the network complexity. As neurons in convolutional layers increase, more subtle features can be captured, which might keep pace with the increasing model output size.

Conclusion
Flexoelectricity is a nonlocal electromechanical coupling phenomenon between strain gradient and electric polarization, which is normally computationally costly for simulating structures with complex topology and further optimizing. We developed an IGA model with the NURBS trimming technique to firstly realize the amiable simulation of porous flexoelectric structures. Then a CNN-based surrogate model is proposed to predict the identical output from our IGA model. Instead of conventional pixel-to-pixel prediction, predicting physical fields of the control points helps to eliminate loss in model accuracy, as well as to save training time of deep learning. Moreover, details of CNN model building are discussed, including the choice of activation functions, dropout layers, and optimizers based on our IGA dataset. Finally, the CNN model performances by changing model output size and convolutional layers size have been analyzed, which might be useful to build similar CNN frameworks predicting nonlocal mech-physical problems.

Appendix 1: Further derivation of IGA formulation
The contribution of a Gaussian integral point to Eq. (11) is where i and W i , are the Jacobian and the weight of the integration point, respectively; and are the corresponding gradient operators for and ; is Hessian operator for . According to the gradient operation of NURBS, , , and can be written as which can all be derived from the explicit expressions of NURBS base functions in Eqs. (8)-(10).

Straight-sided triangular elements
The transformations of the straight-sided triangle element from the physical space to the Gaussian quadrature space is shown in Fig. 18a.

3
The Jacobian of a straight-sided triangular element is determined by transformation R and S , i.e.,

Curved-sided triangular elements
The transformations of the curved-sided triangle element from the physical space to the Gaussian quadrature space is shown in Fig. 18b. The Jacobian of a straight-sided triangular element is determined by transformation P , Q , R and S , i.e., While Transformation S is identical as the one for straight-sided triangular elements, Transformation   Funding Open Access funding enabled and organized by Projekt DEAL.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.