Computational Machine Learning Representation for the Flexoelectricity Effect in Truncated Pyramid Structures

: In this study, machine learning representation is introduced to evaluate the flexoelectricity effect in truncated pyramid nanostructure under compression. A Non-Uniform Rational B-spline (NURBS) based IGA formulation is employed to model the flexoelectricity. We investigate 2D system with an isotropic linear elastic material under plane strain conditions discretized by 45×30 grid of B-spline elements. Six input parameters are selected to construct a deep neural network (DNN) model. They are the Young's modulus, two dielectric permittivity constants, the longitudinal and transversal flexoelectric coefficients and the order of the shape function. The outputs of interest are the strain in the stress direction and the electric potential due flexoelectricity. The dataset are generated from the forward analysis of the flexoelectric model. 80% of the dataset is used for training purpose while the remaining is used for validation by checking the mean squared error. In addition to the input and output layers, the developed DNN model is composed of four hidden layers. The results showed high predictions capabilities of the proposed method with much lower computational time in comparison to the numerical model.


Introduction
In flexoelectric material, the coupling between polarization and strain gradient produces an electro electromechanical coupling due to a mechanical stress or an applied electric field. It exists in several materials. Though, flexoelectricity is possible in different dielectrics, including those with centrosymmetric crystal structures, and is thus a more general electromechanical coupling mechanism than piezoelectricity. Flexoelectricity in solids was introduced first in the 1950s by Mashkevich et al. [Mashkevich and Tolpygo (1957)] but received little attention. Recent developments, however, in nanotechnology have shed the light on flexoelectricity as a size dependent phenomenon due to the large strain gradients that are obtainable at small length scales. The effect of flecxoelectricty is much clear in nanostructures [Nguyen, Mao, Yeh et al. (2013); Krichen and Sharma (2016)]. In this regard, Micro-Nano electromechanical sensors and actuators are increasingly used in numerous applications providing high power density and allowing a broader range of material choice ] and hence micromechanical modeling is essential [Bek, Hamdia, Rabczuk et al. (2018)]. Among the emerging usage of flexoelectricity is the energy harvesters which promise to have one of the widest technological impact replacing the traditional batteries and showing a preferable environmentally alternative [Deng, Kammoun, Erturk et al. (2014)]. In particular, energy harvesters have been used to convert ambient mechanical energy into electrical energy [Priya and Inman (2009)]. Despite the advantages offered by flexoelectricity, research in this field is still in its infancy. The well-known configurations to evaluate the effect of flexoelectricity by generating strain gradients in simple geometries are bending of thin cantilever beam and compression of truncated pyramids. Laboratory tests are limited in measuring, explaining and quantifying the key flexoelectric phenomena at specific conditions [Catalan, Lubk, Vlooswijk et al. (2011);Baskaran, He, Chen et al. (2011);Bhaskar, Banerjee, Abdollahi et al. (2016)]. Besides, further theoretical basis studies have been presented to investigate the flexoelectric phenomena in various structures under different configurations [Sharma, Maranganti and Sharma (2007); Mao and Purohit (2015); Yang, Liang and Shen (2015)]. However, relatively few numerical works have been documented to simulate flexoelectricity. In this regard, fourth order partial differential equations (PDEs) of flexoelectricity has been solved in computational finite element (FE) methods [Abdollahi, Peco, Millán et al. (2014); Abdollahi, Millán, Peco et al. (2015); Mao, Purohit and Aravas (2016)]. A mixed finite element formulation is used to discretize the higher order equations of flexoelectricity [Nanthakumar, Zhuang, Park et al. (2017)]. Besides, mixed FE method was developed using Lagrangian multiplier method to relate the displacement field and its gradient [Deng, Deng, Yu et al. (2017); Deng, Deng and Shen (2018)]. A Non-Uniform Rational B-spline (NURBS) was adopted for the discretization of the solution of fourth order PDEs, where the flexoelectricity effect have been investigated by defining the enhancement in the energy conversion factor [Ghasemi, Park and Rabczuk (2017)]. A distinct approach in this regard is the machine learning representation. Machine learning methods such as artificial neural networks (ANN) has recently proven to be a viable alternative approach to solve complex problems in material design. ANN is stochastic approach based on computational intelligence. It is based on correlating the inputs to the output/s parameters of interest by means of mathematics and statistics methods without the need to an explicit definition of the relationships between inputs and outputs. By using limited amount of training data, ANN has the ability to learn and identify the pattern rapidly. Also, future predictions of the problem can be predicted with much less effort. The classical ANN is composed of three layers: input, one hidden and output layers. Meanwhile, the architecture of deep neural network (DNN) is the same except it has more than two hidden layers. This allows DNN to provide a higher learning representation compared to the classical ANN.
In this study, machine learning approach is presented to gain insight into better understanding the key phenomena governing the flexoelectricity in truncated pyramid structures. DNN is employed for the computational design to relate different input parameters with the predicted behavior aiming at finding low cost constitutive method. Using a set of data for the selected inputs, a numerical model is solved based on the NURBS-based IGA formulation. The input data and the corresponding flexoelectricity material responses form the dataset required to build the model. In the following section, a brief of the material and method used in this study are introduced. Then, Section 3 provides a discussion of the results. Finally, the paper is ended with summary in Section 4.

Materials and method 2.1 Numerical Modeling
The linear continuum theory is well known for modeling flexoelectricity in which a fourth-order flexoelectric coupling tensor is proposed. The electric polarization is a sum of the effects from the dielectric and piezoelectric responses, and the linear polarization response due to the flexoelectricity [Yudin and Tagantsev (2013)]. Mathematically, the relation between the electric polarization, P and the associated strain, , due to mechanical stress is = ⋅ + : + (1) where E is the macroscopic electric field, is the second-order dielectric tensor, is the third order piezoelectric tensor, is the strain tensor, is the fourth-order (including both direct and converse effects) flexoelectric tensor, and ∇ is the spatial gradient of [Sharma, Maranganti and Sharma (2007); Catalan, Sinnamon and Gregg (2004)]. Because of the high-order spatial derivatives term in Eq. (1), conventional finite elements cannot be applied to solve the system. Hence, a Non-Uniform Rational B-spline (NURBS) has been developed recently for discretization in order to solve the continuum equations of flexoelectricity in 2D [Ghasemi, Park and Rabczuk (2017)]. The truncated pyramid under compression load is one of the common configuration in determining the effect of flexoelectricity. We consider two-dimension truncated pyramid geometries with unit width (truncated triangle) by assuming an isotropic linear elasticity under plane strain conditions. Fig. 1 depicts the geometry and the boundary conditions. At the top surface, a uniformly distributed force, F is applied in the negative direction of Y-axis, while the bottom is mechanically fixed. The electric potential is fixed to zero at the top and is constant but unknown at the bottom. The upper face has a length of 1 that is linearly increasing function of the depth up to 2 . Due to the difference in the surface areas of the top and bottom surfaces, the applied force generates different tractions at the top and bottom surfaces, resulting in a strain gradients and accordingly a flexoelectric polarization is generated. A simple analytical formulation to estimate the effective piezoelectric constant accounting for the flexoelectricity in one-dimensional analysis (i.e., only the effect of the longitudinal flexoelectric coefficient, 11 ) was presented by Cross and co-workers [Cross (2006); Zhu, Fu, Li et al. (2006)]. Nevertheless, in the purpose of this study we provide computational modeling based on two dimensional flexoelectricity effect.

Machine learning representation
To model the relation between the inputs and output parameters machine learning modeling can be designed through a predefined training process. The objective of is to build an approximation function ( ) that maps the relation between the inputs ( ) and the outputs ( ). Based on biological learning, ANN in the form of multilayer feedforward networks is adopted in this study. It contains a large number of interconnected processing units called neurons or nodes [Hamdia, Lahmer, Nguyen-Thoi et al. (2015); Hamdia, Arafa and Alqedra (2018)]. These units are grouped together into three or more layers where the neighboring layers are connected by weights forming a large network. The first layer receives information from the input parameters and transmits it to one or several hidden layers, and then evaluates the predictions through the output layer. The network learns by analyzing multiple datasets and adjusting connection weights [Haykin (1999)]. We use neural network with multiple hidden layers called deep neural network (DNN). It uses gradient descent forms for training (updating the weights and bias) via the backward propagation algorithm. Inputs from previous layers are linked to a neuron by the corresponding weights and bias. The weighted sums are applied to an activation function to determine the neuron output, see Fig. 2. The neuron receives data from the preceding layer and consequently proceeds it to the next layer. The tan-sigmoid is considered as the V y x activation function for the hidden layer because it has yielded higher prediction accuracy. The output signal for the -th neuron is expressed as where, is the number of neurons in the preceding layer and −1 is its output signal, while and are the connecting weights and the bias. Finally the signal from the neurons of the last hidden layer is passed to the output layer [Rafiq, Bugmann and Easterbrook (2001)]. The network learns iteratively from several datasets. The predicted outputs are compared with target outputs and accordingly the weights and bias of the neural network are updated to minimize the mean squared error ( ). For training (updating the weights and bias) gradient descent forms are used with Levenberg-Marquardt training algorithm via the error back-propagation (BP) method.

Results
In our analysis, we consider material parameters that fit the behavior of single crystals of barium titanate (BaTiO3). The length of the top and bottom faces are, respectively, 1 =750 μm and 2 =2250 μm discretized by 45×30 grid of B-spline elements. It is made of nonpiezoelectric material obtained through setting 31 =0. The Poisson's ratio ( ) is found to have insignificant effect [Hamdia, Ghasemi, Zhuang et al. (2018)] and hence fixed at 0.37, whereas Young's modulus is taken into account as input parameter. There are also two dielectric permittivity constant 11 and 33 . The longitudinal and transversal flexoelectric coefficients 11 and 12 are considered. Their measurements have shown a high scatter in the experimental investigation which in turn greatly exceed theoretical estimates. We adopt a range of variation of [10-10 5 ] nC/m to define this scatter. In addition to these material parameters, the order of the shape function (p) in the x and y-directions is introduced.  DNN is developed to construct a relation between the addressed input parameters and the corresponding behaviour in each element. The most accurate architecture (minimum mean squared error) is found to be composed of four hidden layers with 30 neurons, 25 neurons, 25 neurons, and 20 neurons, respectively. The DNN model has been constructed by assigning 80% of the data as training data set. The remaining is assigned for testing and validation purpose. Two outputs of are studied using the DNN. They are the distributions of strain in the direction of the applied stress, 22 , and the electric potential, . The computational results and the predictions of the DNN as well are shown in Figure  3 for a randomly selected sample. It illustrates, the distribution of the electric potential due to the mechanical load and the produced strain ( 22 ).
The electric potential and the strain are vary along the entire geometry due to variation of strain gradient that caused by the applied uniform force. Sharp changes of the strain and electric potential near the corners are observed. At the surface subjected to the load (upper face), the electric potential is much lower than those away. The highest values of the electric potential are produced at the inclined surfaces. As expected, also, the upper surface is undergoing compressive strain and decays gradually as going away from it. The predictions of DNN and the computational results are in tune. Nevertheless, the DNN predicts the strain ( 22 ) with higher accuracy. The reason could be that the distribution of the electric potential has higher gradients across the full domain. A more refined discretization for the B-spline elements may improve the results. In spite of the discrepancies in the distribution of the electric potential, the predicted behaviour and results, which have been obtained by a sufficiently fine discretization, are in good agreement with the benchmark example of Abdollahi et al. [Abdollahi, Peco, Millán et al. (2014)] from the point of view both the values and the field distribution.

Summary
In the flexoelectric materials, the polarization is related to not only the strain as in pizoelectric materials, but also to the gradient of strain that enhances their energy conversion efficiency. The truncated pyramid under compression is one of the common configuration to evaluate the flexoelectricity effect in nanostructures. In this set-up, significant strain gradients are generated due to the differences in the areas of the widths of the top and bottom surfaces. In this study, the governing equations of flexoelectricity were modeled by a NURBS-based IGA formulation as ground truth for machine learning representation. We investigated 2D system (truncated triangle) with an isotropic linear elasticity under plane strain conditions. It was discretized by 45×30 grid of B-spline elements. Subsequently, an on-demand property prediction model based deep neural networks (DNN) algorithm was developed to define a mapping between the input parameters and the output of interest. Six input parameters defining the material properties of barium titanate (BaTiO 3) were selected to establish the database. Meanwhile the strain in stress direction and the electric potential due flexoelectricity have been evaluated. The database has been divided randomly into two groups: training and testing datasets. Among variety of architectures that has been trained, a DNN model composed by four hidden layers with 30 neurons, 25 neurons, 25 neurons, and 20 neurons, respectively was found to be the most accurate structure with minimum mean squared error. The predictions of DNN revealed high accuracy predictions compared to the corresponding actual highly demanding computational model. It has been found that at the top surface (with lower area), the electric potential was less than those away and its highest values were produced at the inclined surfaces.