Sensitivity analysis of flexoelectric materials surrogate model based on the isogeometric finite element method

In this paper proposes a sensitivity analysis method based on a Polynomial Chaos Expansion (PCE) surrogate model for flexoelectric materials. The non-uniform rational B-splines (NURBS) basis functions to discretize the fourth-order partial differential equation for flexoelectricity and obtains a deterministic solution (electric potential). The mathematical expressions of surrogate model for the flexoelectric materials are established by considering uncertain parameters such as independent Young’s modulus, concentrated load and flexoelectric constants. The sensitivity expression is found by derivation the mathematical expression for the surrogate model. Moreover the finite difference method (FDM) are conducted in numerical examples to demonstrate the validity and correctness of the proposed algorithm.

1 Introduction

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.

2 The sensitivity expressions of the surrogate model are established by considering three kinds of different material parameters, respectively.

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.

2 Basic formula of polynomial chaos expansion

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

$\begin{array}{c}\hfill f\left(\mathbf{r}\right)=\sum _{\alpha }^{N-1}{\theta }_{\alpha }{\mathrm{\Psi }}_{\alpha }\left(\mathbf{r}\right),\end{array}(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

$\begin{array}{c}\hfill {\mathrm{\Psi }}_{\alpha }\left(\mathbf{r}\right)=\prod _{i=1}^n{{\psi }_{\alpha }}_i\left(r_i\right),\end{array}(2)$

where ${{\psi }_{\alpha }}_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:

$\begin{array}{c}\hfill ⟨{\mathrm{\Psi }}_s\left(\mathbf{r}\right),{\mathrm{\Psi }}_t\left(\mathbf{r}\right)⟩=\int {\mathrm{\Psi }}_s\left(\mathbf{r}\right){\mathrm{\Psi }}_t\left(\mathbf{r}\right)P\left(\mathbf{r}\right)\mathrm{d}\mathbf{r}={\delta }_{st},\end{array}(3)$

where $⟨\cdot ⟩$ 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].

$\mathbf{\Psi }=\left[\begin{array}{ccc}\hfill {\mathrm{\Psi }}_0\left({\mathbf{r}}^1\right)\hfill & \hfill \dots \hfill & \hfill {\mathrm{\Psi }}_{N-1}\left({\mathbf{r}}^1\right)\hfill \\ \hfill \cdot \hfill & \hfill \hfill & \hfill \cdot \hfill \\ \hfill \cdot \hfill & \hfill \hfill & \hfill \cdot \hfill \\ \hfill \cdot \hfill & \hfill \hfill & \hfill \cdot \hfill \\ \hfill {\mathrm{\Psi }}_0\left({\mathbf{r}}^Z\right)\hfill & \hfill \dots \hfill & \hfill {\mathrm{\Psi }}_{N-1}\left({\mathbf{r}}^Z\right)\hfill \end{array}\right]\mathbf{\text{and}}\mathit{\theta }=\left\{\begin{array}{c}\hfill {\theta }_0\hfill \\ \hfill \cdot \hfill \\ \hfill \cdot \hfill \\ \hfill \cdot \hfill \\ \hfill {\theta }_{N-1}\hfill \end{array}\right\},(4)$

where Z is the number of random variables sample points. The PCE coefficient vector can be expressed as

$\mathit{\theta }={\left({\mathbf{\Psi }}^{\mathbf{T}}\mathbf{\Psi }\right)}^{-1}{\mathbf{\Psi }}^{\mathbf{T}}\mathbf{U},(5)$

where $\mathbf{U}={\left\{{\{f(\mathbf{r})\}}_{z=1}^Z\right\}}^{\mathbf{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_{\alpha }\left(r\right)=\frac{1}{{\left(-1\right)}^{\alpha }{e}^{-\frac{r^2}{2}}}\frac{d^{\alpha }}{d{r}^{\alpha }}\left[e^{-\frac{r^2}{2}}\right]=\alpha !\sum _{k=0}^{\left[\alpha /2\right]}{\left(-1\right)}^k\frac{1}{k!2^k\left(\alpha -2k\right)!}{r}^{\alpha -2k},(6)$

where $[\frac{\alpha }{2}]$ is an integer and less than or equal to $\frac{\alpha }{2}$ r denotes the random variable. The three-term recurrence equation for the Hermite orthogonal polynomial is expressed as

$H_{\alpha +1}\left(r\right)=r{H}_{\alpha }\left(r\right)-\alpha H_{\alpha -1}\left(r\right).(7)$

TABLE 1

TABLE 1. The probability distributions of different random variables and corresponding orthogonal polynomials [4].

The first six polynomials of H_α(r) are

$\begin{array}{cc}\hfill & H_{\alpha }\left(r\right)=1\alpha =0,\hfill \\ \hfill & H_{\alpha }\left(r\right)=r\alpha =1,\hfill \\ \hfill & H_{\alpha }\left(r\right)=r^2-1\alpha =2,\hfill \\ \hfill & H_{\alpha }\left(r\right)=r^3-3r\alpha =3,\hfill \\ \hfill & H_{\alpha }\left(r\right)=r^4-6r^2+3\alpha =4,\hfill \\ \hfill & H_{\alpha }\left(r\right)=r^5-10r^3+15r\alpha =5.\hfill \end{array}(8)$

3 IGA discretization of the control equations for flexoelectricity

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

$\begin{array}{ll}\hfill {\int }_{\mathrm{\Omega }}& \left(C_{\mathit{ijkl}}\delta S_{ij}{S}_{kl}-e_{kij}{E}_k\delta S_{ij}-{\mu }_{\mathit{lijk}}{E}_l\delta S_{ij,k}-{\kappa }_{ij}\delta E_i{E}_j-e_{ikl}\delta E_i{S}_{kl}-{\mu }_{ijkl}\delta E_i{S}_{jk,l}\right)d\mathrm{\Omega }-{\int }_{{\mathrm{\Gamma }}_t}{\bar{t}}_i\delta u_{i}d{\mathrm{\Gamma }}_t+{\int }_{{\mathrm{\Gamma }}_D}\omega \delta \phi d{\mathrm{\Gamma }}_D=0,\hfill \end{array}(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; $\overline{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

$\begin{array}{c}\hfill N_{i,0}\left(\xi \right)=\left\{\begin{array}{cc}1\hfill & \mathrm{i}\mathrm{f}{\xi }_i\le \xi <{\xi }_{i+1}\hfill \\ 0\hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\hfill \end{array}\right.,\end{array}(10)$

and for p = 1, 2, 3,

$\begin{array}{c}\hfill N_{i,p}\left(\xi \right)=\frac{\xi -{\xi }_i}{{\xi }_{i+p}-{\xi }_i}{N}_{i,p-1}\left(\xi \right)+\frac{{\xi }_{i+p+1}-\xi }{{\xi }_{i+p+1}-{\xi }_{i+1}}{N}_{i+1,p-1}\left(\xi \right).\end{array}(11)$

Figure 1 illustrates the results of visualizing the B-spline basis functions for both directional knot vectors with ${\mathrm{\Xi }}^1=[\begin{array}{ccccccc}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0.5\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill \end{array}]$ and ${\mathrm{\Xi }}^2=[\begin{array}{ccccccc}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0.5\hfill & \hfill 1\hfill & \hfill 1\hfill & \hfill 1\hfill \end{array}]$. 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

$\begin{array}{c}\hfill \left[\begin{array}{cc}\hfill {\mathbf{A}}_{uu}\hfill & \hfill {\mathbf{A}}_{u\phi }\hfill \\ \hfill {\mathbf{A}}_{\phi u}\hfill & \hfill {\mathbf{A}}_{\phi \phi }\hfill \end{array}\right]\left[\begin{array}{c}\hfill \mathbf{u}\hfill \\ \hfill \mathrm{\Phi }\hfill \end{array}\right]=\left[\begin{array}{c}\hfill {\mathbf{f}}_u\hfill \\ \hfill {\mathbf{f}}_{\phi }\hfill \end{array}\right],\end{array}(12)$

where the matrix corresponding to the displacements is

$\begin{array}{c}\hfill {\mathbf{A}}_{uu}=\sum _e{\int }_{{\mathrm{\Omega }}_e}{\left({\mathbf{B}}_u\right)}^{\mathbf{T}}\mathbf{C}\left({\mathbf{B}}_u\right)\mathrm{d}{\mathrm{\Omega }}_e,\end{array}(13)$

and

$\begin{array}{c}\hfill {\mathbf{B}}_u=\left[\begin{array}{cccccccc}\hfill \frac{\partial N_1}{\partial x}\hfill & \hfill \frac{\partial N_2}{\partial x}\hfill & \hfill \cdots \hfill & \hfill \frac{\partial N_{\mathit{ncp}}}{\partial x}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill & \hfill \frac{\partial N_1}{\partial y}\hfill & \hfill \frac{\partial N_2}{\partial y}\hfill & \hfill \cdots \hfill & \hfill \frac{\partial N_{\mathit{ncp}}}{\partial y}\hfill \\ \hfill \frac{\partial N_1}{\partial y}\hfill & \hfill \frac{\partial N_2}{\partial y}\hfill & \hfill \cdots \hfill & \hfill \frac{\partial N_{\mathit{ncp}}}{\partial y}\hfill & \hfill \frac{\partial N_1}{\partial x}\hfill & \hfill \frac{\partial N_2}{\partial x}\hfill & \hfill \cdots \hfill & \hfill \frac{\partial N_{\mathit{ncp}}}{\partial y}\hfill \end{array}\right],\end{array}(14)$

$\begin{array}{cc}\hfill & \mathbf{C}=\left(\frac{Y}{\left(1+\nu \right)\left(1-2\nu \right)}\right)\left[\begin{array}{ccc}\hfill 1-\nu \hfill & \hfill \nu \hfill & \hfill 0\hfill \\ \hfill \nu \hfill & \hfill 1-\nu \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \left(\frac{1}{2}-\nu \right)\hfill \end{array}\right],\hfill \end{array}(15)$

where ν is Poisson’s ratio and Y is the Young’s modulus. The matrix of displacement and electric field coupling are

$\begin{array}{cc}\hfill & {\mathbf{A}}_{u\phi }=\sum _e{\int }_{{\mathrm{\Omega }}_e}\left[{\left({\mathbf{B}}_u\right)}^{\text{T}}{\mathbf{e}}^{\text{T}}\left({\mathbf{B}}_{\phi }\right)+{\left({\mathbf{H}}_u\right)}^{\text{T}}{\mathit{\mu }}^{\text{T}}\left({\mathbf{B}}_{\phi }\right)\right]\mathrm{d}{\mathrm{\Omega }}_e,\hfill \\ \hfill & {\mathbf{A}}_{\phi u}=\sum _e{\int }_{{\mathrm{\Omega }}_e}\left[{\left({\mathbf{B}}_{\phi }\right)}^{\text{T}}\mathbf{e}\left({\mathbf{B}}_u\right)+{\left({\mathbf{B}}_{\phi }\right)}^{\text{T}}\mathit{\mu }\left({\mathbf{H}}_u\right)\right]\mathrm{d}{\mathrm{\Omega }}_e.\hfill \end{array}(16)$

The electric field and Hessian matrices are

$\begin{array}{cc}\hfill & {\mathbf{B}}_{\phi }=\left[\begin{array}{ccc}\hfill \frac{\partial N_1}{\partial x}\hfill & \hfill \cdots \hfill & \hfill \frac{\partial N_{\mathit{ncp}}}{\partial x}\hfill \\ \hfill \frac{\partial N_1}{\partial y}\hfill & \hfill \cdots \hfill & \hfill \frac{\partial N_{\mathit{ncp}}}{\partial y}\hfill \end{array}\right],\hfill \\ \hfill & {\mathbf{H}}_u=\left[\begin{array}{cccccccc}\hfill \frac{{\partial }^2{N}_1}{\partial x^2}\hfill & \hfill \frac{{\partial }^2{N}_2}{\partial x^2}\hfill & \hfill \cdots \hfill & \hfill \frac{{\partial }^2{N}_{\mathit{ncp}}}{\partial x^2}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill & \hfill \frac{{\partial }^2{N}_1}{\partial y\partial x}\hfill & \hfill \frac{{\partial }^2{N}_2}{\partial y\partial x}\hfill & \hfill \cdots \hfill & \hfill \frac{{\partial }^2{N}_{\mathit{ncp}}}{\partial y\partial x}\hfill \\ \hfill \frac{{\partial }^2{N}_1}{\partial y\partial x}\hfill & \hfill \frac{{\partial }^2{N}_2}{\partial y\partial x}\hfill & \hfill \cdots \hfill & \hfill \frac{{\partial }^2{N}_{\mathit{ncp}}}{\partial y\partial x}\hfill & \hfill \frac{{\partial }^2{N}_1}{\partial x^2}\hfill & \hfill \frac{{\partial }^2{N}_2}{\partial x^2}\hfill & \hfill \cdots \hfill & \hfill \frac{{\partial }^2{N}_{\mathit{ncp}}}{\partial x^2}\hfill \\ \hfill \frac{{\partial }^2{N}_1}{\partial x\partial y}\hfill & \hfill \frac{{\partial }^2{N}_2}{\partial x\partial y}\hfill & \hfill \cdots \hfill & \hfill \frac{{\partial }^2{N}_{\mathit{ncp}}}{\partial x\partial y}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \cdots \hfill & \hfill 0\hfill & \hfill \frac{{\partial }^2{N}_1}{\partial y^2}\hfill & \hfill \frac{{\partial }^2{N}_2}{\partial y^2}\hfill & \hfill \cdots \hfill & \hfill \frac{{\partial }^2{N}_{\mathit{ncp}}}{\partial y^2}\hfill \\ \hfill \frac{{\partial }^2{N}_1}{\partial y^2}\hfill & \hfill \frac{{\partial }^2{N}_2}{\partial y^2}\hfill & \hfill \cdots \hfill & \hfill \frac{{\partial }^2{N}_{\mathit{ncp}}}{\partial y^2}\hfill & \hfill \frac{{\partial }^2{N}_1}{\partial x\partial y}\hfill & \hfill \frac{{\partial }^2{N}_2}{\partial x\partial y}\hfill & \hfill \cdots \hfill & \hfill \frac{{\partial }^2{N}_{\mathit{ncp}}}{\partial x\partial y}\hfill \end{array}\right],\hfill \end{array}(17)$

FIGURE 1

FIGURE 1. The specific B-spline basis functions.

The piezoelectric constants and flexoelectric constants matrices are

$\begin{array}{cc}\hfill & \mathbf{e}=\left[\begin{array}{ccc}\hfill 0\hfill & \hfill 0\hfill & \hfill e_{115}\hfill \\ \hfill e_{311}\hfill & \hfill e_{333}\hfill & \hfill 0\hfill \end{array}\right],\hfill \\ \hfill & \mathit{\mu }=\left[\begin{array}{cccccc}\hfill {\mu }_{11}\hfill & \hfill {\mu }_{12}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill {\mu }_{44}\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill {\mu }_{44}\hfill & \hfill {\mu }_{12}\hfill & \hfill {\mu }_{11}\hfill & \hfill 0\hfill \end{array}\right].\hfill \end{array}(18)$

The matrix corresponding to the electric field is

$\begin{array}{c}\hfill {\mathbf{A}}_{\phi \phi }=-\sum _e{\int }_{{\mathrm{\Omega }}_e}{\left({\mathbf{B}}_{\phi }\right)}^{\mathbf{T}}\mathit{\kappa }\left({\mathbf{B}}_{\phi }\right)\mathrm{d}{\mathrm{\Omega }}_e,\end{array}(19)$

where the permittivity constant matrix is

$\begin{array}{c}\hfill \mathit{\kappa }=\left[\begin{array}{cc}\hfill {\kappa }_{11}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\kappa }_{22}\hfill \end{array}\right].\end{array}(20)$

The force and electrical load vectors are

$\begin{array}{cc}\hfill & {\mathbf{f}}_u=\sum _e{\int }_{{\mathrm{\Gamma }}_{t_e}}{{\mathbf{N}}_u}^{\text{T}}{\mathbf{t}}_{\mathrm{\Gamma }}\mathrm{d}{\mathrm{\Gamma }}_{te},\hfill \\ \hfill & {\mathbf{f}}_{\phi }=-\sum _e{\int }_{{\mathrm{\Gamma }}_{D_e}}{{\mathbf{N}}_{\phi }}^{\text{T}}\omega \mathrm{d}{\mathrm{\Gamma }}_{De},\hfill \end{array}(21)$

The subscript e in Ω_e, Γ_te and Γ_De represents the eth finite element in Eqs. 13, 16, 19, 21.

4 Numercial examples

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.

4.1 Model verification

The cantilever beam model with open-circuit electrical boundary conditions and the top free edge subjected to a concentrated load of 200 μN is depicted in Figure 2. The most commonly used BaTiO₃ material was selected for the cantilever beam model, as summarized in Table 2. The mesh and control point information is demonstrated in Figure 3. The boundary condition of the cantilever beam potential is specified on the right side as 0 V. The FEM uses the traditional Lagrangian basis function, which requires a lot of meshing to achieve higher accuracy, but the processing efficiency is lower. However, IGA-FEM uses NURBS basis functions to reduce the preprocessing time by enhancing the order, which improves the computational efficiency. This is the main motivation for using IGA-FEM in this work.

FIGURE 2

FIGURE 2. Cantilever beam model with open circuit electric boundary conditions.

TABLE 2

TABLE 2. Material parameters for cantilever beam.

FIGURE 3

FIGURE 3. FEM discretization showing the control points as green dots.

The general definition of the electromechanical coupling coefficient ${\mathcal{K}}_{\mathit{eff}}$ is

$\begin{array}{c}\hfill {\mathcal{K}}_{\mathit{eff}}=\sqrt{\frac{{\mathcal{W}}_{\mathit{elec}}}{{\mathcal{W}}_{\mathit{mech}}}}=\sqrt{\frac{{\int }_{\mathrm{\Omega }}{E_i}^{\text{T}}{\kappa }_{ij}E\mathrm{d}\mathrm{\Omega }}{{\int }_{\mathrm{\Omega }}{S_{ij}}^{\text{T}}{C}_{\mathit{ijkl}}{S}_{ij}\mathrm{d}\mathrm{\Omega }}},\end{array}(22)$

For the one-dimensional cantilever beam problem, there are only stresses T₁₁ and electric fields E₂ [31]. The electromechanical coupling coefficient of the cantilever beam problem is [34].

$\begin{array}{c}\hfill {\mathcal{K}}_{\mathit{eff}}=\frac{{\chi }_{33}}{1+{\chi }_{33}}\sqrt{\frac{{\kappa }_{33}}{Y}\left({e_{311}}^2+12{\left(\frac{{\mu }_{12}}{h}\right)}^2\right)},\end{array}(23)$

For comparison purposes, we introduce a normalized expression for the electromechanical coupling coefficient as

$\begin{array}{c}\hfill \bar{\mathcal{K}}=\frac{{\mathcal{K}}_{\mathit{eff}}}{{\mathcal{K}}_{\mathrm{p}\mathrm{i}\mathrm{e}\mathrm{z}\mathrm{o}}},\end{array}(24)$

and

$\begin{array}{c}\hfill {\mathcal{K}}_{\text{piezo}}=\frac{{\chi }_{33}}{1+{\chi }_{33}}\sqrt{\frac{{\kappa }_{33}}{Y}{e_{311}}^2}.\end{array}(25)$

The one-dimensional model can be obtained by setting the Poisson’s ratio, the piezoelectric constant e₃₃₃, flexoelectric constant μ₁₁ in the two-dimensional model to zero. The non-piezoelectric material is obtained by setting e₃₁₁ = 0. The normalized electromechanical coupling coefficients obtained by IGA-FEM and analytical solutions are presented in Tables 3, 4. From the Tables 3, 4, it can be seen that the normalized electromechanical coupling coefficients obtained by IGA-FEM are very close to the analytical solution, and the relative errors are within a small range. The increase in standardized thickness will decrease the mechanical properties (electric potential) of the flexoelectric material.

TABLE 3

TABLE 3. The normalized electromechanical coupling coefficient $\bar{\mathcal{K}}$ of non-piezoelectric materials with changing normalized thickness h′.

TABLE 4

TABLE 4. The normalized electromechanical coupling coefficient $\bar{\mathcal{K}}$ of piezoelectric materials with changing normalized thickness h′.

4.2 PCE surrogate model verification

In this section, the output of the IGA-FEM of the cantilever beam model is used to build the PCE surrogate model. The Young’s modulus, concentrated load and two flexoelectric constants as the random input variables are adopted, respectively. Table 5 lists the mean vlues, coefficients of variation, and ranges of sample points for the different random input variables. In this paper, we have used 500 sample points as input parameters for the random variables of the PCE. The Latin Hypercube Sampling (LHS) random number generation module in Matlab is utilized to obtain sample points of random variables. Owing to the potential distribution at the fixed edge is the most obvious, point A is selected as the reference point when establishing the surrogate model.

TABLE 5

TABLE 5. Definitions and the statistical characteristics of the random input variables.

Figure 4 presents the comparison results of PCE and IGA-FEM for piezoelectric and non-piezoelectric materials with Young’s modulus as random variable. As can be seen from Figure 4, the electric potential decreases as the Young’s modulus increases. The PCE calculation results are basically consistent with the IGA-FEM calculation results, which verifies the effectiveness of the algorithm.

FIGURE 4

FIGURE 4. The electric potential at point A of cantilever beam where the Young’s modulus Y is a random variable.

Figure 5 depicts the electric potential obtained by PCE and IGA-FEM for piezoelectric and non-piezoelectric materials with concentrated load as random variables. As can be seen from Figure 5, the electric potential of point A increases with increasing concentrated load. The electric potential of a piezoelectric material subjected to the same concentrated force F is less than that of a non-piezoelectric material.

FIGURE 5

FIGURE 5. The electric potential at point A of cantilever beam where the concentrated load F is a random variable.

The results of PCE and IGA-FEM for piezoelectric and non-piezoelectric materials are shown in Figure 6, where the random variables are the flexoelectric constants μ₁₁ and μ₁₂. It can be seen from Figure 6 that the electric potential presents reverse changes as the two different flexoelectric constants are changed.

FIGURE 6

FIGURE 6. The electric potential at point A of the cantilever beam for different random variables. (A) The electric potential at point A of the cantilever beam when the random variable is μ₁₁. (B) The electric potential at point A of the cantilever beam when the random variable is μ₁₂.

4.3 Sensitivity analysis

In this section, we perform a sensitivity analysis of the surrogate model for the mechanical properties of flexoelectric materials obtained in Section 4.2. The sensitivity expression of the surrogate model for the mechanical properties of flexoelectric materials without considering piezoelectric effect are

$\begin{array}{cc}\hfill & \frac{d\mathrm{\Phi }\left(Y\right)}{dY}=-2.42\times 10^{-15}+2.48\times 10^{-26}Y-7.35\times 10^{-38}{Y}^2,\hfill \\ \hfill & \frac{d\mathrm{\Phi }\left(F\right)}{dF}=0.36697,\hfill \\ \hfill & \frac{d\mathrm{\Phi }\left({\mu }_{11}\right)}{d{\mu }_{11}}=-39.0240,\hfill \\ \hfill & \frac{d\mathrm{\Phi }\left({\mu }_{12}\right)}{d{\mu }_{12}}=78.4851.\hfill \end{array}(26)$

The sensitivity expression of the surrogate model for the mechanical properties of flexoelectric materials considering piezoelectric effect are

$\begin{array}{cc}\hfill & \frac{d\bar{\mathrm{\Phi }}\left(Y\right)}{dY}=-2.24\times 10^{-15}+2.34\times 10^{-26}Y-6.99\times 10^{-38}{Y}^2,\hfill \\ \hfill & \frac{d\bar{\mathrm{\Phi }}\left(F\right)}{dF}=0.31776,\hfill \\ \hfill & \frac{d\bar{\mathrm{\Phi }}\left({\mu }_{11}\right)}{d{\mu }_{11}}=-36.2533,\hfill \\ \hfill & \frac{d\bar{\mathrm{\Phi }}\left({\mu }_{12}\right)}{d{\mu }_{12}}=81.6950.\hfill \end{array}(27)$

The sensitivity values obtained from the PCE surrogate model are compared with the global finite difference method (FDM) defined by

$\begin{array}{c}\hfill \frac{d\mathrm{\Phi }\left(r\right)}{dr}=\frac{\mathrm{\Phi }\left(r+\mathrm{\Delta }r\right)-\mathrm{\Phi }\left(r\right)}{\mathrm{\Delta }r}.\end{array}(28)$

To investigate the accuracy of DSM and FDM, one gives the relative error of DSM and FDM for different random variable r. The relative error of the sensitivity results obtained by FDM and PCE is

$\begin{array}{c}\hfill {ϵ}_{\mathit{err}}=\frac{\left|{\mathrm{\Phi }}_{\mathit{PCE}}-{\mathrm{\Phi }}_{\mathit{FDM}}\right|}{\left|{\mathrm{\Phi }}_{\mathit{FDM}}\right|}.\end{array}(29)$

Owing to the first sensitivity expressions in Eqs. 27, 28 are a quadratic function of Young’s modulus. Therefore, we need to calculate the sensitivity value of the surrogate model through the mean value of Young’s modulus. The mean value of Young’s modulus is lists in Table 5. The rest of the sensitivity values are directly selected from the derivative results. Tables 6, 7 present a comparison of the sensitivity values of PCE and FDM for non-piezoelectric and piezoelectric materials. From Tables 6, 7, it can be seen that the sensitivity values obtained by PCE and FDM are basically the same, and the relative errors are very small. In terms of the value of sensitivity, the flexoelectric constant μ₁₂ has a greater effect on the electric potential of the flexoelectric material and the material is more sensitive to the flexoelectric constant μ₁₂.

TABLE 6

TABLE 6. Comparison of PCE and FDM of sensitivity results for non-piezoelectric materials.

TABLE 7

TABLE 7. Comparison of PCE and FDM of sensitivity results for non-piezoelectric materials.

5 Conclusion

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.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

Author contributions

Conceptualization, HL; Data curation, HL; Formal analysis, YC and YX; Investigation, XG; Methodology, HL and JZ; Project administration, XY, Software, HL and JZ; Supervision, XY; Validation, YX; Visualization, YC; Writing—original draft, HL and XG. All authors have read and agreed to the published version of the manuscript.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Hamdia KM, Ghasemi H, Zhuang X, Alajlan N, Rabczuk T. Sensitivity and uncertainty analysis for flexoelectric nanostructures. Comput Methods Appl Mech Eng (2018) 337:95–109. doi:10.1016/j.cma.2018.03.016

CrossRef Full Text | Google Scholar

2. Ahmadpoor F, Sharma P. Flexoelectricity in two-dimensional crystalline and biological membranes. Nanoscale (2015) 7:16555–70. doi:10.1039/C5NR04722F

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Yudin PV, Tagantsev AK. Fundamentals of flexoelectricity in solids. Nanotechnology (2013) 24:432001. doi:10.1088/0957-4484/24/43/432001

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Hughes TJ, Cottrell JA, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Comput Methods Appl Mech Eng (2005) 194:4135–95. doi:10.1016/J.CMA.2004.10.008

CrossRef Full Text | Google Scholar

5. Chen L, Wang Z, Peng X, Yang J, Wu P, Lian H. Modeling pressurized fracture propagation with the isogeometric BEM. Geomechanics Geophys Geo-Energy Geo-Resources (2021) 7:51. doi:10.1007/s40948-021-00248-3

CrossRef Full Text | Google Scholar

6. Chen L, Lian H, Xu Y, Li S, Liu Z, Atroshchenko E, et al. Generalized isogeometric boundary element method for uncertainty analysis of time-harmonic wave propagation in infinite domains. Appl Math Model (2023) 114:360–78. doi:10.1016/j.apm.2022.09.030

CrossRef Full Text | Google Scholar

7. Chen L, Lu C, Lian H, Liu Z, Zhao W, Li S, et al. Acoustic topology optimization of sound absorbing materials directly from subdivision surfaces with isogeometric boundary element methods. Comput Methods Appl Mech Eng (2020) 362:112806. doi:10.1016/j.cma.2019.112806

CrossRef Full Text | Google Scholar

8. Chen L, Lian H, Liu Z, Gong Y, Zheng C-J, Bordas S. Bi-material topology optimization for fully coupled structural-acoustic systems with isogeometric FEM-BEM. Eng Anal Boundary Elem (2022) 135:182–95. doi:10.1016/j.enganabound.2021.11.005

CrossRef Full Text | Google Scholar

9. Chen L, Marburg S, Zhao W, Liu C, Chen H. Implementation of isogeometric fast multipole boundary element methods for 2d half-space acoustic scattering problems with absorbing boundary condition. J Theor Comput Acoust (2018) 27:1850024. doi:10.1142/S259172851850024X

CrossRef Full Text | Google Scholar

10. Chen L, Lu C, Zhao W, Chen H, Zheng C-J. Subdivision surfaces—boundary element accelerated by fast multipole for the structural acoustic problem. J Theor Comput Acoust (2020) 28:2050011. doi:10.1142/S2591728520500115

CrossRef Full Text | Google Scholar

11. Chen L, Liu C, Zhao W, Liu L. An isogeometric approach of two dimensional acoustic design sensitivity analysis and topology optimization analysis for absorbing material distribution. Comput Methods Appl Mech Eng (2018) 336:507–32. doi:10.1016/j.cma.2018.03.025

CrossRef Full Text | Google Scholar

12. Chen L, Lian H, Liu Z, Chen H, Atroshchenko E, Bordas S. Structural shape optimization of three dimensional acoustic problems with isogeometric boundary element methods. Comput Methods Appl Mech Eng (2019) 355:926–51. doi:10.1016/j.cma.2019.06.012

CrossRef Full Text | Google Scholar

13. Chen L, Zhang Y, Lian H, Atroshchenko E, Ding C, Bordas SP. Seamless integration of computer-aided geometric modeling and acoustic simulation: Isogeometric boundary element methods based on catmull-clark subdivision surfaces. Adv Eng Softw (2020) 149:102879. doi:10.1016/j.advengsoft.2020.102879

CrossRef Full Text | Google Scholar

14. Chen L, Lian H, Natarajan S, Zhao W, Chen X, Bordas S. Multi-frequency acoustic topology optimization of sound-absorption materials with isogeometric boundary element methods accelerated by frequency-decoupling and model order reduction techniques. Comput Methods Appl Mech Eng (2022) 395:114997. doi:10.1016/j.cma.2022.114997

CrossRef Full Text | Google Scholar

15. Cheng K, Lu Z. Adaptive sparse polynomial chaos expansions for global sensitivity analysis based on support vector regression. Comput Structures (2018) 194:86–96. doi:10.1016/j.compstruc.2017.09.002

CrossRef Full Text | Google Scholar

16. Hurtado J, Barbat A. Monte Carlo techniques in computational stochastic mechanics. Arch Comput Methods Eng (1998) 5:3–29. doi:10.1007/BF02736747

CrossRef Full Text | Google Scholar

17. Chen L, Cheng R, Li S, Lian H, Zheng C, Bordas SP. A sample-efficient deep learning method for multivariate uncertainty qualification of acoustic–vibration interaction problems. Comput Methods Appl Mech Eng (2022) 393:114784. doi:10.1016/j.cma.2022.114784

CrossRef Full Text | Google Scholar

18. Xu Y, Li H, Chen L, Zhao J, Zhang X. Monte Carlo based isogeometric stochastic finite element method for uncertainty quantization in vibration analysis of piezoelectric materials. Mathematics (2022) 10:1840. doi:10.3390/math10111840

CrossRef Full Text | Google Scholar

19. Honda R. Stochastic BEM with spectral approach in elastostatic and elastodynamic problems with geometrical uncertainty. Eng Anal Boundary Elem (2005) 29:415–27. doi:10.1016/j.enganabound.2005.01.007

CrossRef Full Text | Google Scholar

20. Liu WK, Belytschko T, Mani A. Random field finite elements. Int J Numer Methods Eng (1986) 23:1831–45. doi:10.1002/nme.1620231004

CrossRef Full Text | Google Scholar

21. Kamiński M. Stochastic perturbation approach to engineering structure vibrations by the finite difference method. J Sound Vibration (2002) 251:651–70. doi:10.1006/jsvi.2001.3850

CrossRef Full Text | Google Scholar

22. Kamiński M. On generalized stochastic perturbation-based finite element method. Commun Numer Methods Eng (2006) 22:23–31. doi:10.1002/cnm.795

CrossRef Full Text | Google Scholar

23. Zhang B-Y, Ni Y-Q. A hybrid sequential sampling strategy for sparse polynomial chaos expansion based on compressive sampling and Bayesian experimental design. Comput Methods Appl Mech Eng (2021) 386:114130. doi:10.1016/j.cma.2021.114130

CrossRef Full Text | Google Scholar

24. Le Maître O, Knio OM. Spectral methods for uncertainty quantification: With applications to computational fluid dynamics. Berlin, Germany: Springer Science & Business Media (2010).

Google Scholar

25. Blatman G, Sudret B. Sparse polynomial chaos expansions and adaptive stochastic finite elements using a regression approach. Comptes Rendus Mécanique (2008) 336:518–23. doi:10.1016/j.crme.2008.02.013

CrossRef Full Text | Google Scholar

26. Palar PS, Tsuchiya T, Parks GT. Multi-fidelity non-intrusive polynomial chaos based on regression. Comput Methods Appl Mech Eng (2016) 305:579–606. doi:10.1016/j.cma.2016.03.022

CrossRef Full Text | Google Scholar

27. Zhou Y, Lu Z, Cheng K, Shi Y. An expanded sparse Bayesian learning method for polynomial chaos expansion. Mech Syst Signal Process (2019) 128:153–71. doi:10.1016/j.ymssp.2019.03.032

CrossRef Full Text | Google Scholar

28. Saltelli A, Ratto M, Andres T, Campolongo F, Cariboni J, Gatelli D, et al. Global sensitivity analysis: The primer. New York, NY, USA: John Wiley & Sons (2008).

Google Scholar

29. Hamdia KM, Silani M, Zhuang X, He P, Rabczuk T. Stochastic analysis of the fracture toughness of polymeric nanoparticle composites using polynomial chaos expansions. Int J Fracture (2017) 206:215–27. doi:10.1007/s10704-017-0210-6

CrossRef Full Text | Google Scholar

30. Hauseux P, Hale JS, Bordas SP. Accelerating Monte Carlo estimation with derivatives of high-level finite element models. Comput Methods Appl Mech Eng (2017) 318:917–36. doi:10.1016/J.CMA.2017.01.041

CrossRef Full Text | Google Scholar

31. Majdoub MS, Sharma P, Cagin T. Enhanced size-dependent piezoelectricity and elasticity in nanostructures due to the flexoelectric effect. Phys Rev B (2008) 77:125424. doi:10.1103/PhysRevB.77.125424

CrossRef Full Text | Google Scholar

32. Ghasemi H, Park HS, Rabczuk T. A level-set based Iga formulation for topology optimization of flexoelectric materials. Comput Methods Appl Mech Eng (2017) 313:239–58. doi:10.1016/j.cma.2016.09.029

CrossRef Full Text | Google Scholar

33. Ghasemi H, Park HS, Alajlan N, Rabczuk T. A computational framework for design and optimization of flexoelectric materials. Int J Comput Methods (2020) 17:1850097. doi:10.1142/S0219876218500974

CrossRef Full Text | Google Scholar

34. Abdollahi A, Peco C, Millán D, Arroyo M, Arias I. Computational evaluation of the flexoelectric effect in dielectric solids. J Appl Phys (2014) 116:093502. doi:10.1063/1.4893974

CrossRef Full Text | Google Scholar