Flexoelectricity is a new electromechanical energy conversion mechanism that can be an alternative to piezoelectricity [1]. The flexoelectricity is relatively weak in bulk crystalline materials, resulting in little attention. However, with the advancement of nanotechnology, huge strain gradients can be obtained at small-length scales, leading to a new understanding of the flexoelectricity as a size-dependent phenomenon [2]. As compared to piezoelectricity, flexoelectricity theoretically be present in all dielectrics, including those with centrosymmetric crystal structures, and is therefore a more versatile electromechanical coupling mechanism [3]. The traditional Lagrangian interpolation function of Finite Element Method (FEM) can only provide C ⁰ continuity requirements. The C ⁰ continuity of the FEM cannot satisfy the C ¹ continuity requirement of the fourth-order partial differential equation for flexoelectricity. This requires other numerical methods to achieve the C ¹ continuity requirement. Isogeometric analysis (IGA) is one of the most popular numerical methods. It satisfies the continuity of C ¹ by enhancing the order of NURBS basis function [4]. IGA initially developed to unify computer-aided design (CAD) and computer-aided engineering (CAE), but the remarkable characteristics of IGA basis functions such as NURBS has been applied to many applications including mechanics of fracture [5], electromagnetics [6], acoustics [7–10], and optimizations [11–14].

The input parameters of simulation models are often characterized by high uncertainty, and the model parameters are difficult and inaccurate to estimate [15]. This can largely lead us to make erroneous judgments about the issues of concern. The uncertainty quantification of input parameters is an efficient way to address uncertainty, and it examines the uncertainty in the model from the input parameters. Uncertainty analysis methods including Monte Carlo simulation (MCs) [16–18], the random spectral approach [19, 20] and the perturbation technique [21, 22] are frequently used to take into account the impact of uncertainty on the system response. However, with improving accuracy requirements, modeling of target simulation has become extremely complex, and its implementation is costly and time-consuming. The commonly used MCs approach is costly and challenging to implement for uncertainty quantification when many samples and model observations are required [23]. The surrogate modeling approach uses the relationship between inputs and outputs in a basic mathematical model to establish a new method for replacing complex analytical or computational models. The development of surrogate modeling techniques appropriate for solving practical engineering problems provides the required model observations and can reduce the computational cost of uncertainty quantification.

The Polynomial Chaos Expansion (PCE) becomes a prominent alternative modeling method in the field of uncertainty quantification (UQ) with low training cost when modeling extremely complex systems. The main implementation process of PCE is to use several polynomials to expand the response of random variables. The model response is expressed as a polynomial function of the input by determining the PCE coefficients of the polynomial components. These polynomial functions are orthogonal to the probability density functions of the input variables, which makes the calculation easier. The non-intrusive method does not require information about the control equations and is more suitable than the intrusive method for most problems when solving for the PCE coefficients. The non-intrusive methods include projection methods [24] and regression methods [25, 26] in which regression methods are more popular because of their efficiency in dealing with multivariate problems [27]. The sensitivity analysis (SA) quantitatively measures which the uncertainty of different input parameters contributes to the output uncertainty [28]. The sensitivity index is usually used to indicate the influence of each individual input parameter on the output [29]. Some complex problems in practical engineering do not have a definite input-output mathematical expression. It makes it difficult for engineers to perform sensitivity analysis on complex problems. The technique of surrogate modeling, such as PCE, can easily quantify the influence of the input parameters on the output by building mathematical expressions for the mechanical properties of complex problems depending on the inputs and outputs.

This paper lays out a procedure for solving the sensitivity problems of flexoelectric materials. This approach consists of two novel points:

1 The IGA-FEM and PCE are employed to establish a surrogate model for the flexoelectric materials.

The remaining sections of the essay are structured as follows. The foundations of PCE in uncertainty quantification are presented in Section 2. Three introduces the principles of the isogeometric Finite Element Method for the statics of flexoelectric materials. Section 4 validates the IGA-FEM, PCE surrogate models and the sensitivity values of PCE surrogate models of the flexoelectric structure using numerical examples, followed by conclusions in Section 5.

The basic idea of PCE is to replace the system model with an orthogonal polynomial defined by random variables, and then obtain a surrogate model expression by solving for the PCE coefficients. For the system model with n-dimensional independent random variables, the output function of the truncated PCE model with total expand order p can be expressed as f r = ∑ α N − 1 θ α Ψ α r , (1) where the total number N of polynomial terms of order p is computed by N = (n+p)!/(n!p!). θ _α are polynomial coefficients that are unknown. Ψ _α are multivariate orthogonal polynomials defined by the tensor product of univariate orthogonal polynomials as Ψ α r = ∏ i = 1 n ψ α i r i , (2) where ψ α i ( r i ) is a univariate polynomial with respect to the random variable r _i . The inner product of any two functions defined by Ψ _s (r) and Ψ _t (r), and the probability density function P(r) of r is: Ψ s r , Ψ t r = ∫ Ψ s r Ψ t r P r d r = δ s t , (3) where ⋅ is the expectation operator. δ _st is the Kronecker symbol, which is equal to one when s = t and otherwise zero.

In general, the probability density functions for different distributions correspond to different orthogonal polynomials, e.g., Legendre polynomials correspond to uniform distributions and Gaussian distributions correspond to Hermite polynomials. Table 1 lists the common univariate orthogonal polynomials and their corresponding probability distributions. There are various approaches to calculating the polynomial chaos expansion coefficient θ _α . The collocation method and the least-squares minimization problem are commonly employed to obtain the PCE coefficients. The orthogonal matrix and the PCE coefficient vector can be obtained according to Eq. 1 as [30]. Ψ = Ψ 0 r 1 … Ψ N − 1 r 1 ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ Ψ 0 r Z … Ψ N − 1 r Z and θ = θ 0 ⋅ ⋅ ⋅ θ N − 1 , (4) where Z is the number of random variables sample points. The PCE coefficient vector can be expressed as θ = Ψ T Ψ − 1 Ψ T U , (5) where U = { f ( r ) } z = 1 Z T . In this paper, an orthogonal polynomial with random variables satisfying a Gaussian distribution is used as an example. The recurrence relation for the Hermite orthogonal polynomials corresponding to the Gaussian distribution is H α r = 1 − 1 α e − r 2 2 d α d r α e − r 2 2 = α ! ∑ k = 0 α / 2 − 1 k 1 k ! 2 k α − 2 k ! r α − 2 k , (6) where [ α 2 ] is an integer and less than or equal to α 2 r denotes the random variable. The three-term recurrence equation for the Hermite orthogonal polynomial is expressed as H α + 1 r = r H α r − α H α − 1 r . (7)

The first six polynomials of H _α (r) are H α r = 1 α = 0 , H α r = r α = 1 , H α r = r 2 − 1 α = 2 , H α r = r 3 − 3 r α = 3 , H α r = r 4 − 6 r 2 + 3 α = 4 , H α r = r 5 − 10 r 3 + 15 r α = 5 . (8)

In this section, we summarize the controlling equations for dielectric solids considering the flexoelectric effect. More information see [31–33] and references included therein. The weak form of the flexoelectric control equation is ∫ Ω C ijkl δ S i j S k l − e k i j E k δ S i j − μ lijk E l δ S i j , k − κ i j δ E i E j − e i k l δ E i S k l − μ i j k l δ E i S j k , l d Ω − ∫ Γ t t ̄ i δ u i d Γ t + ∫ Γ D ω δ φ d Γ D = 0 , (9) where u _i denotes displacement; φ denotes the electric potential; C _ijkl represents the fourth-order elasticity tensor; the mechanical strain is denoted by S _ij , e _ijk is the third-order piezoelectric tensor, the electric field is defined as E _i = −φ _i ; the fourth-order total flexoelectric tensor is denoted by μ _ijkl ; κ _ij is dielectric tensor of second order; t i ̄ denotes the mechanical traction; ω is surface charge density. The physical domain is denoted by Ω, with boundaries Γ _t and Γ _D corresponding to mechanical traction and electric displacements, respectively. In order to obtain the FEM form of the governing Eq. 9, the B-spline basis function N _i,p (ξ) obtained by recursively defining the Cox-de-Boor formula is N i , 0 ξ = 1 i f ξ i ≤ ξ < ξ i + 1 0 o t h e r w i s e , (10) and for p = 1, 2, 3, N i , p ξ = ξ − ξ i ξ i + p − ξ i N i , p − 1 ξ + ξ i + p + 1 − ξ ξ i + p + 1 − ξ i + 1 N i + 1 , p − 1 ξ . (11)

Figure 1 illustrates the results of visualizing the B-spline basis functions for both directional knot vectors with Ξ 1 = [ 0 0 0 0.5 1 1 1 ] and Ξ 2 = [ 0 0 0 0.5 1 1 1 ] . The N _i,p and N _i,q are 2^nd order basis functions. The richness of B-spline basis function can be intuitively seen in Figure 1. The richness of the basis functions provides the groundwork for solving the fourth-order partial differential equation for flexoelectricity. Using Eqs. 10, 11 to discretize Eq. 9, the linear algebraic system of equations for the flexoelectricity control equation is obtained as A u u A u φ A φ u A φ φ u Φ = f u f φ , (12) where the matrix corresponding to the displacements is A u u = ∑ e ∫ Ω e B u T C B u d Ω e , (13) and B u = ∂ N 1 ∂ x ∂ N 2 ∂ x ⋯ ∂ N ncp ∂ x 0 0 0 ⋯ 0 0 ⋯ 0 ∂ N 1 ∂ y ∂ N 2 ∂ y ⋯ ∂ N ncp ∂ y ∂ N 1 ∂ y ∂ N 2 ∂ y ⋯ ∂ N ncp ∂ y ∂ N 1 ∂ x ∂ N 2 ∂ x ⋯ ∂ N ncp ∂ y , (14) C = Y 1 + ν 1 − 2 ν 1 − ν ν 0 ν 1 − ν 0 0 0 1 2 − ν , (15) where ν is Poisson’s ratio and Y is the Young’s modulus. The matrix of displacement and electric field coupling are A u φ = ∑ e ∫ Ω e B u T e T B φ + H u T μ T B φ d Ω e , A φ u = ∑ e ∫ Ω e B φ T e B u + B φ T μ H u d Ω e . (16) The electric field and Hessian matrices are B φ = ∂ N 1 ∂ x ⋯ ∂ N ncp ∂ x ∂ N 1 ∂ y ⋯ ∂ N ncp ∂ y , H u = ∂ 2 N 1 ∂ x 2 ∂ 2 N 2 ∂ x 2 ⋯ ∂ 2 N ncp ∂ x 2 0 0 ⋯ 0 0 0 ⋯ 0 ∂ 2 N 1 ∂ y ∂ x ∂ 2 N 2 ∂ y ∂ x ⋯ ∂ 2 N ncp ∂ y ∂ x ∂ 2 N 1 ∂ y ∂ x ∂ 2 N 2 ∂ y ∂ x ⋯ ∂ 2 N ncp ∂ y ∂ x ∂ 2 N 1 ∂ x 2 ∂ 2 N 2 ∂ x 2 ⋯ ∂ 2 N ncp ∂ x 2 ∂ 2 N 1 ∂ x ∂ y ∂ 2 N 2 ∂ x ∂ y ⋯ ∂ 2 N ncp ∂ x ∂ y 0 0 ⋯ 0 0 0 ⋯ 0 ∂ 2 N 1 ∂ y 2 ∂ 2 N 2 ∂ y 2 ⋯ ∂ 2 N ncp ∂ y 2 ∂ 2 N 1 ∂ y 2 ∂ 2 N 2 ∂ y 2 ⋯ ∂ 2 N ncp ∂ y 2 ∂ 2 N 1 ∂ x ∂ y ∂ 2 N 2 ∂ x ∂ y ⋯ ∂ 2 N ncp ∂ x ∂ y , (17)

The piezoelectric constants and flexoelectric constants matrices are e = 0 0 e 115 e 311 e 333 0 , μ = μ 11 μ 12 0 0 0 μ 44 0 0 μ 44 μ 12 μ 11 0 . (18)

The matrix corresponding to the electric field is A φ φ = − ∑ e ∫ Ω e B φ T κ B φ d Ω e , (19) where the permittivity constant matrix is κ = κ 11 0 0 κ 22 . (20)

The force and electrical load vectors are f u = ∑ e ∫ Γ t e N u T t Γ d Γ t e , f φ = − ∑ e ∫ Γ D e N φ T ω d Γ D e , (21) The subscript e in Ω _e , Γ _te and Γ _De represents the eth finite element in Eqs. 13, 16, 19, 21.

In this section, we verify the accuracy of the IGA-FEM for solving the fourth-order partial differential equation using a benchmark example of a cantilever beam. After that, some random variable sample points are selected to obtain the output of IGA-FEM, which is used to build a polynomial chaos expansion surrogate model. Finally, the sensitivity results of the surrogate model for the mechanical properties of flexoelectric materials are verified by several numerical examples. For the cantilever beam model, we postulate that the model satisfies plane strain linear elastic isotropy.

In this paper, a sensitivity analysis method of surrogate model based on isogeometric stochastic Finite Element Method is proposed for flexoelectric materials. The NURBS basis functions with high-order continuity are used to discretize the fourth-order partial differential equation for flexoelectricity. The Polynomial Chaos Expansion (PCE) is utilized to develop a surrogate model for the mechanical properties of flexoelectric materials. The sensitivity values of the surrogate model are obtained by considering three kinds of different parameters, respectively. Numerical examples illustrate the flexoelectric material is more sensitive to the flexoelectric constant μ ₁₂. Additionally, the current technology will also be used to three-dimensional piezoelectric and flexoelectric problems.