A New Meshless "Fragile Points Method (FPM)" Based on A Galerkin Weak-Form for 2D Flexoelectric Analysis
1 A NEW MESHLESS โFRAGILE POINTS METHOD (FPM)โ BASED ON A GALERKIN WEAK-FORM FOR 2D FLEXOELECTRIC ANALYSIS
Yue Guan , Leiting Dong , Satya N. Atluri Texas Tech University, Lubbock, TX Beihang University, Beijing, China
ABSTRACT
A meshless Fragile Points Method (FPM) is presented for analyzing 2D flexoelectric problems. Local, simple, polynomial and discontinuous trial and test functions are generated with the help of a local meshless differential quadrature approximation of the first three derivatives. Interior Penalty Numerical Fluxes are employed to ensure the consistency of the method. Based on a Galerkin weak-form formulation, the present FPM leads to symmetric and sparse matrices, and avoids the difficulties of numerical integration in the previous meshfree methods. Numerical examples including isotropic and anisotropic materials with flexoelectric and piezoelectric effects are provided as validations. The present method is much simpler than the Finite Element Method, or the Element-Free Galerkin (EFG) and Meshless Local Petrov-Galerkin (MLPG) methods, and the numerical integration of the weak form is trivially simple.
Keywords: Fragile Points Method (FPM), meshless method, flexoelectricity, gradient elasticity
1. INTRODUCTION
The flexoelectricity describes an electromechanical coupling effect between the electric polarization and mechanical strain gradient [1-3]. As a result of the recent miniaturized trend of electromechanical systems, which leads a tremendous growth in the scale of strain gradient, there has been an increasing demand for reliable and accurate theories and numerical methods for flexoelectric analysis [2,4]. In 2006, Maranganti et al. [5] proposed the first continuum theory for flexoelectricity, in which the internal energy density is given as a function of strain, strain gradient, electric polarization and polarization gradient. After that, a number of advanced continuum theories have been carried out, considering the Maxwell stress [6], the nonlinear flexoelectric behavior under larger deformation [7], etc. Generally, the flexoelectric behavior is governed by a fourth-order partial differential equation. As a result, the main difficulty in modelling flexoelectricity for classic element-based methods lies in the ๐ถ๐ถ continuity requirement. On one hand, the Finite Element Method (FEM) based on Argyris triangular element [8] can be employed to remedy the ๐ถ๐ถ requirement. Sladek et al. [9] have also proposed a ๐ถ๐ถ continuous FEM analyzing the effects of electric field and strain gradients. However, both of the approaches require 9 degrees of freedom (DOFs) at each node, and thus make it computationally expensive. Alternatively, some mixed FEM [3,4] using Lagrangian multipliers as independent variables are also available, requiring only ๐ถ๐ถ continuity. Yet they have even more DOFs in each element, making it prohibitively expensive for practical use. On the other hand, meshless methods like the Local Maximum-Entropy (LME) Meshfree Method [10], Element-Free Galerkin (EFG) Method [11], and Meshless Local Petrov-Galerkin (MLPG) Method [12] have also shown their capability in analyzing flexoelectric or strain-gradient effects. In the EFG and MLPG method, with the help of Moving Least Square (MLS) approximation, global shape and trial functions with ๐ถ๐ถ continuity can be generated. However, these trial functions based on MLS are rational and considerably complicated. As a result, their numerical integration in the weak form in either EFG or MLPG can be tedious. Actually, the difficulty in domain integration is a common problem in many Galerkin meshfree methods [13]. The simplest choice is direct nodal integration. Though efficient, and no background mesh required, the direct nodal integration may have stability issues and can be less accurate or non-convergent. A number of modified nodal integration methods are proposed, in order to ensure the accuracy and stability, which in turn sacrifice the efficiency [13]. Some other newer-type quadrature schemes are also adopted [14]. Alternatively, the difficulty of numerical integration can be solved fundamentally by employing simple shape functions like polynomials. The Fragile Points Method (FPM), in contrast to most of the previous meshless methods, is based on simple, polynomial, piecewise-continuous trial and test functions. Consequently, the classic Gaussian quadrature is applied, and numerical integrals are relatively simple, especially when simple subdomain shapes are used (e.g., quadrilateral or triangle). Numerical Flux Corrections are employed to ensure the continuity condition. On the other hand, when compared with the element-based methods, the FPM is generated using point stiffness matrices. Thus, it has benefits in inserting or removing points and bypassing the mesh distortion problem. The discontinuity of the trial function also presents a great potential in analyzing systems involving cracks, ruptures and fragmentations. Furthermore, when compared with the EFG and MLPG methods, the FPM also has advantages like symmetric matrices, and the shape functions pass through the nodal data directly. The FPM has already shown great accuracy and efficiency in solving 2D thermal conductivity and elasticity problems. In this study, we formulate and apply the FPM for analyzing linear flexoelectric problems. 2
2. FLEXOELECTRICITY THEORY AND BOUNDARY-VALUE PROBLEM
Phenomenologically, the flexoelectricity describes an electric polarization generated by the mechanical strain-gradient: ๐๐ ๐๐ = ๐๐ฬ ๐๐๐๐๐๐๐๐ ๐ ๐ ๐๐๐๐๐๐ , (1) where ๐๐ = [ ๐๐ , ๐๐ ] T is the electrical polarization vector. ๐ฎ๐ฎ = [ ๐ข๐ข , ๐ข๐ข ] T is the displacement vector in a 2D domain ๐บ๐บ . ๐ ๐ ๐๐๐๐๐๐ = ๐ข๐ข ๐๐ , ๐๐๐๐ is the second gradient of displacement (i.e., the strain-gradient), and ๐๐ฬ ๐๐๐๐๐๐๐๐ are elements of a fourth order flexoelectric tensor ( ๐๐ , ๐๐ , ๐๐ , ๐๐ = 1,2 ). Einstein summation convention is used here. In the current study, the linear flexoelectric theory proposed by [5] is employed, in which the internal energy is described as a function of strain ๐๐ ๐๐๐๐ , strain gradient ๐ ๐ ๐๐๐๐๐๐ and polarization ๐๐ ๐๐ , i.e. ๐๐๏ฟฝ = ๐๐๏ฟฝ ( ๐๐ , ๐๐ , ๐๐ ) . The corresponding governing equations can be written as: ๏ฟฝ๐๐ ๐๐๐๐ โ ๐๐ ๐๐๐๐๐๐ , ๐๐ ๏ฟฝ , ๐๐ + ๐๐ ๐๐ = 0, (2) โ๐๐ ๐๐ , ๐๐๐๐ + ๐๐ ๐๐ , ๐๐ = ๐๐ (3) where ๐๐ is the Cauchy stress, ๐๐ is the double-stress (conjugate of the strain gradient ๐ฟ๐ฟ ), ๐๐ is the electrical potential, ๐๐ is the body force per volume, ๐๐ is the free charge per volume, and ๐๐ is the permittivity of free space. The corresponding boundary conditions are ๐ข๐ข ๐๐ = ๐ข๐ข๏ฟฝ ๐๐ , ๐๐๐๐ ๐๐๐บ๐บ ๐ข๐ข , (4) ๏ฟฝ๐๐ ๐๐๐๐ โ ๐๐ ๐๐๐๐๐๐ , ๐๐ ๏ฟฝ๐๐ ๐๐ + ๏ฟฝ๏ฟฝโ ๐๐๐ก๐ก ๐๐ ๐๐ ๏ฟฝ๐๐ ๐๐ โ โ ๐๐๐ก๐ก ๏ฟฝ ( ๐๐ ๐๐ ๐๐ ๐๐๐๐๐๐ ) = ๐๐๏ฟฝ ๐๐ , ๐๐๐๐ ๐๐๐บ๐บ ๐๐ , (5) โ ๐๐ ๐ข๐ข ๐๐ = ๐๐ฬ ๐๐ , ๐๐๐๐ ๐๐๐บ๐บ ๐๐ , (6) ๐๐ ๐๐ ๐๐ ๐๐ ๐๐ ๐๐๐๐๐๐ = ๐ ๐ ๏ฟฝ ๐๐ , ๐๐๐๐ ๐๐๐บ๐บ ๐ ๐ , (7) ๐๐ = ๐๐๏ฟฝ , ๐๐๐๐ ๐๐๐บ๐บ ๐๐ , (8) ๏ฟฝโ๐๐ ๐๐ , ๐๐ + ๐๐ ๐๐ ๏ฟฝ๐๐ ๐๐ = โ๐๐๏ฟฝ , ๐๐๐๐ ๐๐๐บ๐บ ๐๐ , (9) where ๏ฟฝ๐ฎ๐ฎ๏ฟฝ , ๐๐๏ฟฝ , ๐๐ฬ , ๐๐๏ฟฝ , ๐๐๏ฟฝ , ๐๐๏ฟฝ๏ฟฝ are known functions, โ ๐๐ = ๐๐ โ ๐๐ is the normal derivative, ๐๐ ๐ก๐ก = ๐๐ โ ๐๐โ ๐๐ the โsurface gradientโ on ๐๐๐บ๐บ . ๐๐๐บ๐บ ๐ข๐ข โช ๐๐๐บ๐บ ๐๐ = ๐๐๐บ๐บ ๐๐ โช ๐๐๐บ๐บ ๐ ๐ = ๐๐๐บ๐บ ๐๐ โช ๐๐๐บ๐บ ๐๐ = ๐๐๐บ๐บ , and ๐๐๐บ๐บ ๐ข๐ข โฉ ๐๐๐บ๐บ ๐๐ = ๐๐๐บ๐บ ๐๐ โฉ ๐๐๐บ๐บ ๐ ๐ = ๐๐๐บ๐บ ๐๐ โฉ ๐๐๐บ๐บ ๐๐ = โ . We consider a material with both piezoelectric and flexoelectric effects. The constitutive equations are: ๐๐ = ๐๐ ฯฮต ๐๐ โ ๐๐ ๐๐ โ ๐๐๐๐ , (10) ๐๐ = ๐๐ ฮผฮบ ๐๐ โ ๐๐ ๐๐ โ ๐๐ ๐๐ (11) ๐๐ = ( ๐๐๏ฟฝ โ ๐๐ ๐๐ ) ๐๐ + ๐๐ T ๐๐ + ๐๐ ๐๐ , (12) where ๐๐ = [ ๐๐ , ๐๐ , 2 ๐๐ ] T , ๐๐ = [ ๐๐ , ๐๐ , ๐๐ ] T , ๐๐ = [ ๐ ๐ , ๐ ๐ , 2 ๐ ๐ , 2 ๐ ๐ , ๐ ๐ , ๐ ๐ ] T , ๐๐ = [ ๐๐ , ๐๐ , ๐๐ , ๐๐ , ๐๐ , ๐๐ ] T , ๐๐ = [ ๐ธ๐ธ , ๐ธ๐ธ ] T , and ๐๐ ๐๐๐๐ = 12 ๏ฟฝ๐ข๐ข ๐๐ , ๐๐ + ๐ข๐ข ๐๐ , ๐๐ ๏ฟฝ , ๐๐ ๐๐๐๐ = ๐๐ ๐๐๐๐ , ๐ ๐ ๐๐๐๐๐๐ = ๐ ๐ ๐๐๐๐๐๐ , ๐๐ ๐๐๐๐๐๐ = ๐๐ ๐๐๐๐๐๐ , ๐ธ๐ธ ๐๐ = โ๐๐ , ๐๐ . The constitutive matrices: ๐๐ ฯฮต = ๐๐๏ฟฝ ฯฮต โ ๐๐ ( ๐๐๏ฟฝ โ ๐๐ ๐๐ ) โ1 ๐๐ T , (13) ๐๐ ฮผฮบ = ๐๐๏ฟฝ ฮผฮบ โ ๐๐ ( ๐๐๏ฟฝ โ ๐๐ ๐๐ ) โ1 ๐๐ , (14) ๐๐ = ๐๐ ( ๐๐๏ฟฝ โ ๐๐ ๐๐ ) โ1 ๐๐ . (15) ๐๐๏ฟฝ ฯฮต , ๐๐๏ฟฝ ฮผฮบ , ๐๐๏ฟฝ , ๐๐ and ๐๐ are matrices of material properties (see Appendix A). ๐๐ is a unit matrix. In the numerical calculations, for simplicity, we only concentrate on materials with isotropic elasticity, cubic symmetry for the flexoelectric tensor, and tetragonal symmetry for the piezoelectric tensor.
3. TRIAL AND TEST FUNCTIONS 3.1 Points and domain partition
First, a set of random points are scattered in the domain. The entire domain can then be partitioned into several nonoverlapping subdomains. Within each subdomain, only one point exists. Multiple partitions are possible here, e.g., the Voronoi Diagram partition (see Fig. 1(a)), quadrilateral partition (see Fig. 1(b)), and triangular partition. The traditional FEM meshing can also be employed. The elements can be converted into subdomains, while the internal point is defined as the center of mass of each subdomain. Thus, the preprocessing module of commercial FEA software like ABAQUS can be helpful for generating the points and subdomains in this proposed FPM. However, unlike the FEM, the trial and test functions in the present FPM are point-based . In order to describe the high-order behavior, the trial function ๐ฎ๐ฎ โ in each subdomain can be written in the form of a third-order Taylor expansion at the corresponding internal point. For instance, in subdomain ๐ธ๐ธ with internal point ๐๐ : ๐ฎ๐ฎ โ ( ๐ฅ๐ฅ , ๐ฆ๐ฆ ) = ๐ฎ๐ฎ + ( ๐ฅ๐ฅ โ ๐ฅ๐ฅ ) ๐ฎ๐ฎ , ๏ฟฝ ๐๐ + ( ๐ฆ๐ฆ โ ๐ฆ๐ฆ ) ๐ฎ๐ฎ , ๏ฟฝ ๐๐ + ( ๐ฅ๐ฅ โ ๐ฅ๐ฅ ) ๐ฎ๐ฎ , ๏ฟฝ ๐๐ + ( ๐ฅ๐ฅ โ ๐ฅ๐ฅ )( ๐ฆ๐ฆ โ ๐ฆ๐ฆ ) ๐ฎ๐ฎ , ๏ฟฝ ๐๐ + ( ๐ฆ๐ฆ โ ๐ฆ๐ฆ ) ๐ฎ๐ฎ , ๏ฟฝ ๐๐ + ( ๐ฅ๐ฅ โ ๐ฅ๐ฅ ) ๐ฎ๐ฎ , ๏ฟฝ ๐๐ + ( ๐ฅ๐ฅ โ ๐ฅ๐ฅ ) ( ๐ฆ๐ฆ โ ๐ฆ๐ฆ ) ๐ฎ๐ฎ , ๏ฟฝ ๐๐ + ( ๐ฅ๐ฅ โ ๐ฅ๐ฅ )( ๐ฆ๐ฆ โ๐ฆ๐ฆ ) ๐ฎ๐ฎ , ๏ฟฝ ๐๐ + ( ๐ฆ๐ฆ โ ๐ฆ๐ฆ ) ๐ฎ๐ฎ , ๏ฟฝ ๐๐ , (16) where ( ๐ฅ๐ฅ , ๐ฆ๐ฆ ) are the coordinates of the point ๐๐ . ๐ฎ๐ฎ is the value of ๐ฎ๐ฎ โ at ๐๐ . The first three derivatives [ ๐๐๏ฟฝ๐ฎ๐ฎ ]| ๐๐ = ๏ฟฝ๐ฎ๐ฎ , , ๐ฎ๐ฎ , , ๐ฎ๐ฎ , , ๐ฎ๐ฎ , , ๐ฎ๐ฎ , , ๐ฎ๐ฎ , , ๐ฎ๐ฎ , , ๐ฎ๐ฎ , , ๐ฎ๐ฎ , ๏ฟฝ T ๏ฟฝ ๐๐ are yet unknown. Here we introduce the local differential quadrature method to determine the high-order derivatives. (A) (B) FIGURE 1:
THE DOMAIN ฮฉ AND ITS PARTITIONS. (A) VORONOI DIAGRAM PARTITION. (B) QUADRILATERAL PARTITION.
The conventional differential quadrature (DQ) method is a widely used numerical discretization technique [15]. However, usually based on 1D basis functions (e.g., polynomials), the conventional DQ method can only be applied along a mesh-line and cannot be incorporated directly with meshless schemes. In recent years, several multi-dimensional DQ methods have been developed based on the meshless Radial Basis Functions (RBFs). In the current work, we employ the local RBF-DQ method proposed by Shu et al. [15] in 2003 which approximates the derivatives directly using a limited number of supporting points. In the current methodology, as shown in Fig. 1, for a given point ๐๐ โ ๐ธ๐ธ (red), its support is defined to involve all its nearest (blue) and second (green) neighboring points. Here the nearest neighboring points are the points in subdomains sharing boundaries with ๐ธ๐ธ , while the second neighboring points are nearest neighbors of the former. For the points on or close to the domain boundary ๐๐๐บ๐บ , the third neighboring points are also taken into consideration to remedy the lack of effective supporting points. These supporting points are named as ๐๐ , ๐๐ , โฆ , ๐๐ ๐๐ . Many RBFs can be used as basis function in the local RBF-DQ method. The multiquadric (MQ), inverse-MQ and Gaussian RBFs are all available in conjunction with the FPM and achieve similar accuracy with appropriate parameters. Here we use the MQ-RBF as an example: ๐๐ ( ๐๐ ) = ๏ฟฝ๐๐ + ๐๐ , (17) where ๐๐ is the radial from the conference point. ๐๐ = ๐๐ ๐ท๐ท is a constant parameter, where ๐ท๐ท is the diameter of the minimal circle enclosing all the supporting points. Thus, any partial derivative of displacement at point ๐๐ can be approximated by a weighted linear sum of the value of ๐ฎ๐ฎ โ at all its supporting points: ๐๐ ๐ ๐ +๐ก๐ก ๐๐๐ฅ๐ฅ ๐ ๐ ๐๐๐ฆ๐ฆ ๐ก๐ก ๐ฎ๐ฎ โ ๏ฟฝ ๐๐ = ๏ฟฝ ๐๐ ๐๐ ( ๐ ๐ , ๐ก๐ก ) ๐ฎ๐ฎ โ | ๐๐ ๐๐ ๐๐๐๐=0 , (18) where ๐๐ ๐๐ ( ๐ ๐ , ๐ก๐ก ) is the weighting coefficient corresponding to point ๐๐ ๐๐ . After rearrangement, the derivatives under study at point ๐๐ can be approximated as: [ ๐๐๏ฟฝ ๐ฎ๐ฎ ]| ๐๐ = ( ๐๐๏ฟฝโจ๐๐ ร ) ๐ฎ๐ฎ E , (19) where ๐ฎ๐ฎ E = [ ๐ข๐ข , ๐ข๐ข , ๐ข๐ข , ๐ข๐ข , โฆ , ๐ข๐ข , ๐ข๐ข ] T . (20) ๐ฎ๐ฎ ๐๐ = ๏ฟฝ ๐ข๐ข ๐๐ , ๐ข๐ข ๐๐ ๏ฟฝ T is the value of ๐ฎ๐ฎ โ at ๐๐ ๐๐ . ๐๐๏ฟฝ is the matrix of weighting coefficients ๐๐ ๐๐ ( ๐ ๐ , ๐ก๐ก ) , ๐๐ ร is a unit matrix, and โจ donates the Kronecker product. The weighting coefficient matrix ๐๐๏ฟฝ can be obtained from:
๐๐๏ฟฝ = ๐๐ โ1 [ ๐๐๏ฟฝ ๐๐ ], (21) where ๐๐ = ๏ฟฝ ๐๐ ( ๐ฅ๐ฅ , ๐ฆ๐ฆ ) ๐๐ ( ๐ฅ๐ฅ , ๐ฆ๐ฆ ) โฆ ๐๐ ( ๐ฅ๐ฅ ๐๐ , ๐ฆ๐ฆ ๐๐ ) โฎ โฎ โฑ โฎ๐๐ ๐๐ ( ๐ฅ๐ฅ , ๐ฆ๐ฆ ) ๐๐ ๐๐ ( ๐ฅ๐ฅ , ๐ฆ๐ฆ ) โฆ ๐๐ ๐๐ ( ๐ฅ๐ฅ ๐๐ , ๐ฆ๐ฆ ๐๐ ) ๏ฟฝ , [ ๐๐๏ฟฝ ๐๐ ] = โฃโขโขโขโขโขโขโขโขโขโก ๐๐ , ( ๐ฅ๐ฅ , ๐ฆ๐ฆ ) โฆ ๐๐ ๐๐ , ( ๐ฅ๐ฅ , ๐ฆ๐ฆ )0 ๐๐ , ( ๐ฅ๐ฅ , ๐ฆ๐ฆ ) โฆ ๐๐ ๐๐ , ( ๐ฅ๐ฅ , ๐ฆ๐ฆ )0 ๐๐ , ( ๐ฅ๐ฅ , ๐ฆ๐ฆ ) โฆ ๐๐ ๐๐ , ( ๐ฅ๐ฅ , ๐ฆ๐ฆ )0 ๐๐ , ( ๐ฅ๐ฅ , ๐ฆ๐ฆ ) โฆ ๐๐ ๐๐ , ( ๐ฅ๐ฅ , ๐ฆ๐ฆ )0 ๐๐ , ( ๐ฅ๐ฅ , ๐ฆ๐ฆ ) โฆ ๐๐ ๐๐ , ( ๐ฅ๐ฅ , ๐ฆ๐ฆ )0 ๐๐ , ( ๐ฅ๐ฅ , ๐ฆ๐ฆ ) โฆ ๐๐ ๐๐ , ( ๐ฅ๐ฅ , ๐ฆ๐ฆ )0 ๐๐ , ( ๐ฅ๐ฅ , ๐ฆ๐ฆ ) โฆ ๐๐ ๐๐ , ( ๐ฅ๐ฅ , ๐ฆ๐ฆ )0 ๐๐ , ( ๐ฅ๐ฅ , ๐ฆ๐ฆ ) โฆ ๐๐ ๐๐ , ( ๐ฅ๐ฅ , ๐ฆ๐ฆ )0 ๐๐ , ( ๐ฅ๐ฅ , ๐ฆ๐ฆ ) โฆ ๐๐ ๐๐ , ( ๐ฅ๐ฅ , ๐ฆ๐ฆ ) โฆโฅโฅโฅโฅโฅโฅโฅโฅโฅโค T , ๐๐ ๐๐ ( ๐ฅ๐ฅ , ๐ฆ๐ฆ ) = ๏ฟฝ ( ๐ฅ๐ฅ โ ๐ฅ๐ฅ ๐๐ ) + ( ๐ฆ๐ฆ โ ๐ฆ๐ฆ ๐๐ ) + ๐๐ โ ๏ฟฝ ( ๐ฅ๐ฅ โ ๐ฅ๐ฅ ) + ( ๐ฆ๐ฆ โ ๐ฆ๐ฆ ) + ๐๐ . Note that the accuracy of local RBF-DQ approximation is excellent for the first derivative. Whereas for higher-order approximations, the accuracy decreases yet are still acceptable. This implies that a mixed formulation may help to improve the accuracy of the current method.
Substitute the approximation of derivatives into Eqn. (16). The trial function ๐ฎ๐ฎ โ in subdomain ๐ธ๐ธ is achieved: ๐ฎ๐ฎ โ ( ๐ฅ๐ฅ , ๐ฆ๐ฆ ) = ๐ต๐ต ( ๐ฅ๐ฅ , ๐ฆ๐ฆ ) ๐ฎ๐ฎ E , (22) where ๐๐ ( ๐ฅ๐ฅ , ๐ฆ๐ฆ ) = ๏ฟฝ๐๐๏ฟฝ โ ๐๐๏ฟฝ + [
๐๐ ๐๐ โฆ ] ร( ๐๐+1 ) ๏ฟฝ โจ๐๐ ร , ๐๐๏ฟฝ = ๏ฟฝ๐ฅ๐ฅ โ ๐ฅ๐ฅ , ๐ฆ๐ฆ โ ๐ฆ๐ฆ , ( ๐ฅ๐ฅ โ ๐ฅ๐ฅ ) , ( ๐ฅ๐ฅ โ ๐ฅ๐ฅ )( ๐ฆ๐ฆ โ๐ฆ๐ฆ ), ( ๐ฆ๐ฆ โ ๐ฆ๐ฆ ) , ( ๐ฅ๐ฅ โ ๐ฅ๐ฅ ) , ( ๐ฅ๐ฅ โ ๐ฅ๐ฅ ) ( ๐ฆ๐ฆ โ๐ฆ๐ฆ ), ( ๐ฅ๐ฅ โ ๐ฅ๐ฅ )( ๐ฆ๐ฆ โ ๐ฆ๐ฆ ) , ( ๐ฆ๐ฆ โ ๐ฆ๐ฆ ) ๏ฟฝ . The trial function and shape function are defined following the same process in each subdomain. Since no continuity requirement is employed at the internal boundaries, the shape and trial functions can be discontinuous. Figure 2(A) shows the graph of a shape function in a unit domain with 50 random points and Voronoi Diagram partition. The trial function simulating an exponential function ๐ข๐ข ๐๐ = ๐๐ โ10๏ฟฝ ( ๐ฅ๐ฅโ0 . ) + ( ๐ฆ๐ฆโ0 . ) is presented in Fig. 2(B). As can be seen, the trial function is a cubic polynomial in each subdomain. It is piecewise-continuous. However, as the shape function in each subdomain also depends on the values of the neighboring points, the entire shape function shows a weakly continuous tendency. The shape and trial functions for the electric potential ๐๐ โ are generated similarly. The test function for displacement ๐ฏ๐ฏ โ and electrical potential ๐๐ โ in the Galerkin weak-form in the FPM are prescribed to possess the same shape as ๐ฎ๐ฎ โ and ๐๐ โ respectively. As a result of the discontinuous trial functions, the FPM also has a great potential in analyzing systems involving cracks, ruptures and fragmentations. (A) (B) FIGURE 2:
THE SHAPE FUNCTION AND TRIAL FUNCTION.
4. Weak-form formulation and Numerical Flux Corrections
We multiply the governing equations (2) and (3) by the test functions ๐ฃ๐ฃ ๐๐โ and ๐๐ ๐๐ respectively in each subdomain ๐ธ๐ธ , then apply the Gauss divergence theorem: ๏ฟฝ ๐ฃ๐ฃ ๐๐ , ๐๐ ๐๐ ๐๐๐๐ ๐ธ๐ธ d ๐บ๐บ + ๏ฟฝ ๐ฃ๐ฃ ๐๐ , ๐๐๐๐ ๐๐ ๐๐๐๐๐๐ ๐ธ๐ธ d ๐บ๐บ = ๏ฟฝ ๐ฃ๐ฃ ๐๐ ๐ธ๐ธ ๐๐ ๐๐ d ๐บ๐บ + ๏ฟฝ ๐ฃ๐ฃ ๐๐ ๏ฟฝ๐๐ ๐๐๐๐ โ ๐๐ ๐๐๐๐๐๐ , ๐๐ ๏ฟฝ๐๐ ๐๐ ๐๐๐ธ๐ธ d ๐ค๐ค + ๏ฟฝ ๐ฃ๐ฃ ๐๐ , ๐๐ ๐๐ ๐๐๐๐๐๐ ๐๐ ๐๐ ๐๐๐ธ๐ธ d ๐ค๐ค (23) ๏ฟฝ ๐๐ ๐๐ , ๐๐ ๐๐ , ๐๐ ๐ธ๐ธ d ๐บ๐บ โ ๏ฟฝ ๐๐ , ๐๐ ๐๐ ๐๐ ๐ธ๐ธ d ๐บ๐บ = ๏ฟฝ ๐๐๐๐ ๐ธ๐ธ d ๐บ๐บ โ ๏ฟฝ ๐๐๏ฟฝโ๐๐ ๐๐ , ๐๐ + ๐๐ ๐๐ ๏ฟฝ๐๐ ๐๐ ๐๐๐ธ๐ธ d ๐ค๐ค (24) where ๐๐๐ธ๐ธ is the boundary of the subdomain, ๐ง๐ง = [ ๐๐ , ๐๐ ] T is the unit vector outward to ๐๐๐ธ๐ธ . Let ๐ค๐ค โ donates the set of all internal boundaries, and ๐ค๐ค = ๐ค๐ค โ + ๐๐๐บ๐บ = ๐ค๐ค โ + ๐๐๐บ๐บ ๐ข๐ข + ๐๐๐บ๐บ ๐๐ = ๐ค๐ค โ + ๐๐๐บ๐บ ๐๐ + ๐๐๐บ๐บ ๐ ๐ = ๐ค๐ค โ + ๐๐๐บ๐บ ๐๐ + ๐๐๐บ๐บ ๐๐ . By summing the equations over the entire domain, we rewrite them in the matrix forms: ๏ฟฝ ๐๐ ๐๐ ( ๐ฏ๐ฏ ) ๐๐ ( ๐ฎ๐ฎ , ๐๐ ) ๐บ๐บ d ๐บ๐บ + ๏ฟฝ ๐๐ ๐๐ ( ๐ฏ๐ฏ ) ๐๐ ( ๐ฎ๐ฎ , ๐๐ ) ๐บ๐บ d ๐บ๐บ โ ๏ฟฝ โฆ๐๐๏ฟฝ T ( ๐ฏ๐ฏ ) โง๐ง๐ง๏ฟฝ { ๐๐ ( ๐ฎ๐ฎ , ๐๐ )} ๐ค๐ค โ d ๐ค๐ค โ ๏ฟฝ { ๐๐๏ฟฝ T ( ๐ฏ๐ฏ )} ๐ง๐ง๏ฟฝ โฆ๐๐ ( ๐ฎ๐ฎ , ๐๐ ) โง ๐ค๐ค โ d ๐ค๐ค โ ๏ฟฝ โฆ๐ฏ๐ฏ ๐๐ โง ๏ฟฝ๐ง๐ง๏ฟฝ ๐๐ ( ๐ฎ๐ฎ , ๐๐ ) โ ๐ง๐ง๏ฟฝ ๐๐ , ( ๐ฎ๐ฎ , ๐๐ ) โ๐ง๐ง๏ฟฝ ๐๐ , ( ๐ฎ๐ฎ , ๐๐ ) ๏ฟฝ d ๐ค๐ค ๐ค๐ค โ โ ๏ฟฝ { ๐ฏ๐ฏ ๐๐ } ๏ฟฝ๐ง๐ง๏ฟฝ ๐๐ ( ๐ฎ๐ฎ , ๐๐ ) โ ๐ง๐ง๏ฟฝ ๐๐ , ( ๐ฎ๐ฎ , ๐๐ ) โ๐ง๐ง๏ฟฝ ๐๐ , ( ๐ฎ๐ฎ , ๐๐ ) ๏ฟฝ ๐ค๐ค โ d ๐ค๐ค = ๏ฟฝ ๐ฏ๐ฏ ๐๐ ๐๐ ๐บ๐บ d ๐บ๐บ + ๏ฟฝ โฆ๐๐๏ฟฝ T ( ๐ฏ๐ฏ ) โง๐ง๐ง๏ฟฝ { ๐๐ ( ๐ฎ๐ฎ , ๐๐ )} โฮฉ d ๐ค๐ค + ๏ฟฝ โฆ๐ฏ๐ฏ ๐๐ โง ๏ฟฝ๐ง๐ง๏ฟฝ ๐๐ ( ๐ฎ๐ฎ , ๐๐ ) โ ๐ง๐ง๏ฟฝ ๐๐ , ( ๐ฎ๐ฎ , ๐๐ ) โ๐ง๐ง๏ฟฝ ๐๐ , ( ๐ฎ๐ฎ , ๐๐ ) ๏ฟฝ ๐๐๐บ๐บ d ๐ค๐ค (25) ๏ฟฝ ๐๐ ๐๐ ๐๐ ( ๐๐ ) ๐๐ ( ๐๐ ) ๐บ๐บ d ๐บ๐บ + ๏ฟฝ ๐๐ ๐๐ ( ๐๐ ) ๐๐ ( ๐ฎ๐ฎ , ๐๐ ) ๐บ๐บ d ๐บ๐บ + ๏ฟฝ โฆ๐๐โง๐ง๐ง ๐๐๐๐ { ๐๐ ๐๐ ( ๐๐ ) + ๐๐ ( ๐ฎ๐ฎ , ๐๐ )} ๐ค๐ค โ d ๐ค๐ค + ๏ฟฝ { ๐๐ } ๐ง๐ง ๐๐๐๐ โฆ๐๐ ๐๐ ( ๐๐ ) + ๐๐ ( ๐ฎ๐ฎ , ๐๐ ) โง ๐ค๐ค โ d ๐ค๐ค = ๏ฟฝ ๐๐๐๐ ๐บ๐บ d ๐บ๐บ โ ๏ฟฝ โฆ๐๐โง๐ง๐ง ๐๐ ๏ฟฝ ๐๐ ๐๐ ( ๐๐ )+ ๐๐ ( ๐ฎ๐ฎ , ๐๐ ) ๏ฟฝ ๐๐๐บ๐บ d ๐ค๐ค (26) where ๐๐๏ฟฝ ( ๐ฎ๐ฎ ) = ๏ฟฝ๐ข๐ข , , ๐ข๐ข , , ๐ข๐ข , , ๐ข๐ข , ๏ฟฝ T . Each internal boundary ๐๐ โ ๐ค๐ค โ is shared by two subdomains, i.e., ๐๐ = ๐๐๐ธ๐ธ โฉ ๐๐๐ธ๐ธ . ๐ง๐ง ๐๐ is defined as the unit vector normal to ๐๐ and pointing outward from ๐ธ๐ธ , that is, ๐ง๐ง ๐๐ = ๐ง๐ง | ๐๐๐ธ๐ธ = โ๐ง๐ง | ๐๐๐ธ๐ธ . The jump operator โฆโโง and average operator { โ } are defined as (for โ๐ค๐ค ): โฆ๐ค๐คโง = ๏ฟฝ ๐ค๐ค | ๐๐๐ธ๐ธ โ ๐ค๐ค | ๐๐๐ธ๐ธ ๐๐ โ ๐ค๐ค โ ๐ค๐ค | ๐๐ ๐๐ โ ๐๐ฮฉ , { ๐ค๐ค } = ๏ฟฝ
12 ( ๐ค๐ค | ๐๐๐ธ๐ธ + ๐ค๐ค | ๐๐๐ธ๐ธ ) ๐๐ โ ๐ค๐ค โ ๐ค๐ค | ๐๐ ๐๐ โ ๐๐ฮฉ . The following boundary term in Eqn. (25) can be integrated by parts: ๏ฟฝ โฆ๐๐๏ฟฝ T ( ๐ฏ๐ฏ ) โง๐ง๐ง๏ฟฝ { ๐๐ ( ๐ฎ๐ฎ , ๐๐ )} โฮฉ d ๐ค๐ค = ๏ฟฝ ๐๐๏ฟฝ T ( ๐ฏ๐ฏ ) ๐ง๐ง๏ฟฝ ๐ง๐ง๏ฟฝ ๐ง๐ง๏ฟฝ ๐๐ ( ๐ฎ๐ฎ , ๐๐ ) ๐๐๐บ๐บ ๐๐ โช๐๐๐บ๐บ ๐ ๐ d ๐ค๐ค โ ๏ฟฝ ๐ฏ๐ฏ T ๐๐ฬ ๐ฌ๐ฌฬ ๐ฌ๐ฌฬ ๐ง๐ง๏ฟฝ ๐๐ , ( ๐ฎ๐ฎ , ๐๐ ) ๐๐๐บ๐บ ๐ข๐ข โช๐๐๐บ๐บ ๐๐ d ๐ค๐ค โ ๏ฟฝ ๐ฏ๐ฏ T ๐๐ฬ ๐ฌ๐ฌฬ ๐ฌ๐ฌฬ ๐ง๐ง๏ฟฝ ๐๐ , ( ๐ฎ๐ฎ , ๐๐ ) ๐๐๐บ๐บ ๐ข๐ข โช๐๐๐บ๐บ ๐๐ d ๐ค๐ค . (27) The corresponding matrices in Eqn. (25) โ (27) are shown in Appendix B. When ( ๐ฎ๐ฎ , ๐๐ ) are exact solutions, due to continuity conditions, โฆ๐๐ ( ๐ฎ๐ฎ , ๐๐ ) โง = , ๏ฟฝ๐ง๐ง๏ฟฝ ๐๐ ( ๐ฎ๐ฎ , ๐๐ ) โ ๐ง๐ง๏ฟฝ ๐๐ , ( ๐ฎ๐ฎ , ๐๐ ) โ๐ง๐ง๏ฟฝ ๐๐ , ( ๐ฎ๐ฎ , ๐๐ ) ๏ฟฝ = and โฆ๐๐ ๐๐ ( ๐๐ ) + ๐๐ ( ๐ฎ๐ฎ , ๐๐ ) โง = for โ๐๐ โ๐ค๐ค โ . Hence, we can eliminate the corresponding terms in Eqn. (25) and (26). Similarly, โฆ๐ฎ๐ฎโง = , โฆ๐๐๏ฟฝ ( ๐ฎ๐ฎ ) โง = and โฆ๐๐โง = 0 . By adding antithetic terms, the previous equations can be written in a symmetric form. Furthermore, in order to impose the essential boundary conditions ( ๐๐๐บ๐บ ๐ข๐ข , ๐๐๐บ๐บ ๐๐ and ๐๐๐บ๐บ ๐๐ ) and the continuity condition across subdomains, we employ the Numerical Flux Correction, which is widely used in Discontinuous Galerkin (DG) Methods [16] to ensure the consistency and stability of the method. Here we use the Interior Penalty (IP) Numerical Flux Corrections. A set of penalty parameters are introduced, in which ๐๐ =[ ๐๐ , ๐๐ , ๐๐ ] and ๐๐ = [ ๐๐ , ๐๐ , ๐๐ ] are related to the essential boundary conditions and the interior continuity conditions respectively. The FPM is only stable when the penalty parameters are large enough. Finally, after substituting the constitutive equations and the following boundary conditions (see Eqn. (4) โ (9)): ๏ฟฝ ๐ฏ๐ฏ ๐๐ ๏ฟฝ ๐ง๐ง๏ฟฝ ๐๐ ( ๐ฎ๐ฎ , ๐๐ ) โ๏ฟฝ๐ง๐ง๏ฟฝ ๐๐ฬ ๐ฌ๐ฌฬ ๐ฌ๐ฌฬ ๐ง๐ง๏ฟฝ ๏ฟฝ๐๐ , ( ๐ฎ๐ฎ , ๐๐ ) โ๏ฟฝ๐ง๐ง๏ฟฝ โ ๐๐ฬ ๐ฌ๐ฌฬ ๐ฌ๐ฌฬ ๐ง๐ง๏ฟฝ ๏ฟฝ๐๐ , ( ๐ฎ๐ฎ , ๐๐ ) ๏ฟฝ ๐๐๐บ๐บ ๐๐ d ๐ค๐ค = ๏ฟฝ ๐ฏ๐ฏ ๐๐ ๐๐๏ฟฝ ๐๐๐บ๐บ ๐๐ d ๐ค๐ค , (28) ๏ฟฝ ๐๐๏ฟฝ T ( ๐ฏ๐ฏ ) ๐ง๐ง๏ฟฝ ๐ง๐ง๏ฟฝ ๐ง๐ง๏ฟฝ ๐๐ ( ๐ฎ๐ฎ , ๐๐ )d ๐ค๐ค ๐๐๐บ๐บ ๐ ๐ = ๏ฟฝ ๐๐๏ฟฝ T ( ๐ฏ๐ฏ ) ๐ง๐ง๏ฟฝ ๐๐๏ฟฝ d ๐ค๐ค ๐๐๐บ๐บ ๐ ๐ , (29) ๏ฟฝ โฆ๐๐โง๐ง๐ง ๐๐ { ๐๐ ๐๐ ( ๐๐ ) + ๐๐ ( ๐ฎ๐ฎ , ๐๐ )} ๐๐๐บ๐บ ๐๐ d ๐ค๐ค = โ ๏ฟฝ ๐๐๐๐๏ฟฝ ๐๐๐บ๐บ ๐๐ d ๐ค๐ค , (30) we can obtain the final symmetric formula for the FPM for flexoelectric problems: ๏ฟฝ ๐๐ T ( ๐ฏ๐ฏ ) ๐๐ ฯฮต ๐๐ ( ๐ฎ๐ฎ ) ๐บ๐บ d ๐บ๐บ + ๏ฟฝ ๐๐ T ( ๐ฏ๐ฏ ) ๐๐ ฮผฮบ ๐๐ ( ๐ฎ๐ฎ ) ๐บ๐บ d ๐บ๐บ โ ๏ฟฝ ๏ฟฝ๐๐ T ( ๐ฏ๐ฏ ) ๐๐ ๐๐ ( ๐ฎ๐ฎ ) + ๐๐ T ( ๐ฏ๐ฏ ) ๐๐ ๐๐ ( ๐ฎ๐ฎ ) ๏ฟฝ ๐บ๐บ d ๐บ๐บ โ ๏ฟฝ ๏ฟฝ๐๐ T ( ๐ฏ๐ฏ ) ๐๐ + ๐๐ T ( ๐ฏ๐ฏ ) ๐๐ ๏ฟฝ๐๐ ( ๐๐ ) ๐บ๐บ d ๐บ๐บ + ๏ฟฝ โโโโโ โฆ๐ฏ๐ฏ T โง ๏ฟฝ ๐ง๐ง๏ฟฝ ๐๐ โ๐ง๐ง๏ฟฝ ๐๐ ๐๐ฬ โ๐ง๐ง๏ฟฝ ๐๐ ๐๐ฬ ๏ฟฝ { ๐๐ ( ๐ฎ๐ฎ )}+{ ๐๐ T ( ๐ฏ๐ฏ )} ๏ฟฝ ๐๐ ๐ง๐ง๏ฟฝ โ๐๐ฬ ๐๐ ๐ง๐ง๏ฟฝ โ๐๐ฬ ๐๐ ๐ง๐ง๏ฟฝ ๏ฟฝ โฆ๐ฎ๐ฎโงโ โโโโ d ๐ค๐ค ๐ค๐ค โ โช๐๐๐บ๐บ ๐ข๐ข โ ๏ฟฝ ๏ฟฝ โฆ๐ฏ๐ฏ T โง๐ง๐ง๏ฟฝ ๐๐ ๏ฟฝ๐๐ , ( ๐๐ ) ๏ฟฝ + โฆ๐ฏ๐ฏ T โง๐ง๐ง๏ฟฝ ๐๐ ๏ฟฝ๐๐ , ( ๐๐ ) ๏ฟฝ๏ฟฝ d ๐ค๐ค ๐ค๐ค โ โช๐๐๐บ๐บ ๐ข๐ข โ ๏ฟฝ ๏ฟฝ โฆ๐ฏ๐ฏ T โง๐ง๐ง๏ฟฝ ๐๐ ฯฮต { ๐๐ ( ๐ฎ๐ฎ )}+{ ๐๐ T ( ๐ฏ๐ฏ )} ๐๐ ฯฮต ๐ง๐ง๏ฟฝ โฆ๐ฎ๐ฎโง๏ฟฝ d ๐ค๐ค ๐ค๐ค โ โช๐๐๐บ๐บ ๐ข๐ข + ๏ฟฝ ๏ฟฝ โฆ๐ฏ๐ฏ T โง๐ง๐ง๏ฟฝ ๐๐ ฮผฮบ ๏ฟฝ๐๐ , ( ๐ฎ๐ฎ ) ๏ฟฝ + ๏ฟฝ๐๐ , ( ๐ฏ๐ฏ ) ๏ฟฝ๐๐ ฮผฮบ ๐ง๐ง๏ฟฝ โฆ๐ฎ๐ฎโง๏ฟฝ d ๐ค๐ค ๐ค๐ค โ โช๐๐๐บ๐บ ๐ข๐ข + ๏ฟฝ ๏ฟฝ โฆ๐ฏ๐ฏ T โง๐ง๐ง๏ฟฝ ๐๐ ฮผฮบ ๏ฟฝ๐๐ , ( ๐ฎ๐ฎ ) ๏ฟฝ + ๏ฟฝ๐๐ , ( ๐ฏ๐ฏ ) ๏ฟฝ๐๐ ฮผฮบ ๐ง๐ง๏ฟฝ โฆ๐ฎ๐ฎโง๏ฟฝ d ๐ค๐ค ๐ค๐ค โ โช๐๐๐บ๐บ ๐ข๐ข + ๏ฟฝ โฆ๐ฏ๐ฏ T โง๐ง๐ง๏ฟฝ ๐๐ { ๐๐ ( ๐๐ )}d ๐ค๐ค ๐ค๐ค โ โช๐๐๐บ๐บ ๐ข๐ข โ ๏ฟฝ ๏ฟฝ โฆ๐๐๏ฟฝ T ( ๐ฏ๐ฏ ) โง๐ง๐ง๏ฟฝ ๐๐ ฮผฮบ { ๐๐ ( ๐ฎ๐ฎ )}+{ ๐๐ T ( ๐ฏ๐ฏ )} ๐๐ ฮผฮบ ๐ง๐ง๏ฟฝ โฆ๐๐๏ฟฝ ( ๐ฎ๐ฎ ) โง๏ฟฝ ๐ค๐ค โ d ๐ค๐ค + ๏ฟฝ ๏ฟฝ โฆ๐๐๏ฟฝ T ( ๐ฏ๐ฏ ) โง๐ง๐ง๏ฟฝ ๐๐ { ๐๐ ( ๐ฎ๐ฎ )}+{ ๐๐ T ( ๐ฏ๐ฏ )} ๐๐ ๐ง๐ง๏ฟฝ โฆ๐๐๏ฟฝ ( ๐ฎ๐ฎ ) โง๏ฟฝ ๐ค๐ค โ d ๐ค๐ค + ๏ฟฝ โฆ๐๐๏ฟฝ T ( ๐ฏ๐ฏ ) โง๐ง๐ง๏ฟฝ ๐๐ { ๐๐ ( ๐๐ )} ๐ค๐ค โ d ๐ค๐ค + ๏ฟฝ ๏ฟฝ ๐ฏ๐ฏ T ๐๐ฬ ๐ฌ๐ฌฬ ๐ฌ๐ฌฬ ๐ง๐ง๏ฟฝ ๐๐ ฮผฮบ ๐๐ , ( ๐ฎ๐ฎ )+ ๐๐ , ( ๐ฏ๐ฏ ) ๐๐ ฮผฮบ ๐ง๐ง๏ฟฝ ๐ฌ๐ฌฬ ๐ฌ๐ฌฬ ๐๐ฬ ๐ฎ๐ฎ๏ฟฝ ๐๐๐บ๐บ ๐ข๐ข d ๐ค๐ค + ๏ฟฝ ๏ฟฝ ๐ฏ๐ฏ T ๐๐ฬ ๐ฌ๐ฌฬ ๐ฌ๐ฌฬ ๐ง๐ง๏ฟฝ ๐๐ ฮผฮบ ๐๐ , ( ๐ฎ๐ฎ )+ ๐๐ , ( ๐ฏ๐ฏ ) ๐๐ ฮผฮบ ๐ง๐ง๏ฟฝ ๐ฌ๐ฌฬ ๐ฌ๐ฌฬ ๐๐ฬ ๐ฎ๐ฎ๏ฟฝ ๐๐๐บ๐บ ๐ข๐ข d ๐ค๐ค โ ๏ฟฝ โโโ ๐ฏ๐ฏ T ๏ฟฝ ๐๐ฬ ๐ฌ๐ฌฬ ๐ฌ๐ฌฬ ๐ง๐ง๏ฟฝ ๐๐ ๐๐ฬ + ๐๐ฬ ๐ฌ๐ฌฬ ๐ฌ๐ฌฬ ๐ง๐ง๏ฟฝ ๐๐ ๐๐ฬ ๏ฟฝ ๐๐ ( ๐ฎ๐ฎ )+ ๐๐ T ( ๐ฏ๐ฏ ) ๏ฟฝ ๐๐ฬ ๐๐ ๐ง๐ง๏ฟฝ ๐ฌ๐ฌฬ ๐ฌ๐ฌฬ ๐๐ฬ + ๐๐ฬ ๐๐ ๐ง๐ง๏ฟฝ ๐ฌ๐ฌฬ ๐ฌ๐ฌฬ ๐๐ฬ ๏ฟฝ ๐ฎ๐ฎโ โโ ๐๐๐บ๐บ ๐ข๐ข d ๐ค๐ค โ ๏ฟฝ ๏ฟฝ ๐ฏ๐ฏ T ๐๐ฬ ๐ฌ๐ฌฬ ๐ฌ๐ฌฬ ๐ง๐ง๏ฟฝ ๐๐ ๐๐ , ( ๐๐ )+ ๐ฏ๐ฏ T ๐๐ฬ ๐ฌ๐ฌฬ ๐ฌ๐ฌฬ ๐ง๐ง๏ฟฝ ๐๐ ๐๐ , ( ๐๐ ) ๏ฟฝ ๐๐๐บ๐บ ๐ข๐ข d ๐ค๐ค โ ๏ฟฝ ๏ฟฝ ๐๐๏ฟฝ T ( ๐ฏ๐ฏ ) ๐ง๐ง๏ฟฝ ๐ง๐ง๏ฟฝ ๐ง๐ง๏ฟฝ ๐๐ ฮผฮบ ๐๐ ( ๐ฎ๐ฎ )+ ๐๐ T ( ๐ฏ๐ฏ ) ๐๐ ฮผฮบ ๐ง๐ง๏ฟฝ ๐ง๐ง๏ฟฝ ๐ง๐ง๏ฟฝ ๐๐๏ฟฝ ( ๐ฎ๐ฎ ) ๏ฟฝ ๐๐๐บ๐บ ๐๐ d ๐ค๐ค + ๏ฟฝ ๏ฟฝ ๐๐๏ฟฝ T ( ๐ฏ๐ฏ ) ๐ง๐ง๏ฟฝ ๐ง๐ง๏ฟฝ ๐ง๐ง๏ฟฝ ๐๐ ๐๐ ( ๐ฎ๐ฎ )+ ๐๐ T ( ๐ฏ๐ฏ ) ๐๐ ๐ง๐ง๏ฟฝ ๐ง๐ง๏ฟฝ ๐ง๐ง๏ฟฝ ๐๐๏ฟฝ ( ๐ฎ๐ฎ ) ๏ฟฝ ๐๐๐บ๐บ ๐๐ d ๐ค๐ค + ๏ฟฝ ๐๐๏ฟฝ T ( ๐ฏ๐ฏ ) ๐ง๐ง๏ฟฝ ๐ง๐ง๏ฟฝ ๐ง๐ง๏ฟฝ ๐๐ ๐๐ ( ๐๐ ) ๐๐๐บ๐บ ๐๐ d ๐ค๐ค โ ๏ฟฝ ๏ฟฝ { ๐๐ T ( ๐ฏ๐ฏ )} ๐๐ +{ ๐๐ T ( ๐ฏ๐ฏ )} ๐๐ ๏ฟฝ ๐ง๐ง ๐๐ โฆ๐๐โง d ๐ค๐ค ๐ค๐ค โ โช๐๐๐บ๐บ ๐๐ + ๏ฟฝ ๐๐ โ ๐๐ ๐๐๐บ๐บ ๐ข๐ข ๐ฏ๐ฏ T ๐ฎ๐ฎ d ๐ค๐ค + ๏ฟฝ ๐๐ โ ๐๐ ๐ค๐ค โ โฆ๐ฏ๐ฏ T โงโฆ๐ฎ๐ฎโง d ๐ค๐ค + ๏ฟฝ ๐๐ ๐๐๐บ๐บ ๐๐ โ ๐๐ ๐๐๏ฟฝ T ( ๐ฏ๐ฏ ) ๐ง๐ง๏ฟฝ ๐ง๐ง๏ฟฝ ๐๐๏ฟฝ ( ๐ฎ๐ฎ )d ๐ค๐ค + ๏ฟฝ ๐๐ ๐ค๐ค โ โ ๐๐ โฆ๐๐๏ฟฝ T ( ๐ฏ๐ฏ ) โง๐ง๐ง๏ฟฝ ๐ง๐ง๏ฟฝ โฆ๐๐๏ฟฝ ( ๐ฎ๐ฎ ) โง d ๐ค๐ค = ๏ฟฝ ๐ฏ๐ฏ T ๐บ๐บ ๐๐ d ๐บ๐บ + ๏ฟฝ ๐ฏ๐ฏ T ๐๐๏ฟฝ ๐๐๐บ๐บ ๐๐ d ๐ค๐ค + ๏ฟฝ ๐๐๏ฟฝ T ( ๐ฏ๐ฏ ) ๐ง๐ง๏ฟฝ ๐๐๏ฟฝ d ๐ค๐ค ๐๐๐บ๐บ ๐ ๐ โ ๏ฟฝ ๐๐ T ( ๐ฏ๐ฏ ) ๐๐ ฯฮต ๐ง๐ง๏ฟฝ ๐ฎ๐ฎ๏ฟฝ d ๐ค๐ค ๐๐๐บ๐บ ๐ข๐ข + ๏ฟฝ ๐๐ T ( ๐ฏ๐ฏ ) ๏ฟฝ๐๐ ๐ง๐ง๏ฟฝ โ ๐๐ฬ ๐๐ ๐ง๐ง๏ฟฝ โ๐๐ฬ ๐๐ ๐ง๐ง๏ฟฝ ๏ฟฝ ๐ฎ๐ฎ๏ฟฝ d ๐ค๐ค ๐๐๐บ๐บ ๐ข๐ข + ๏ฟฝ ๏ฟฝ ๐๐ , ( ๐ฏ๐ฏ ) ๐๐ ฮผฮบ ๏ฟฝ๐ง๐ง๏ฟฝ + ๐ง๐ง๏ฟฝ ๐ฌ๐ฌฬ ๐ฌ๐ฌฬ ๐๐ฬ ๏ฟฝ + ๐๐ , ( ๐ฏ๐ฏ ) ๐๐ ฮผฮบ ๏ฟฝ๐ง๐ง๏ฟฝ + ๐ง๐ง๏ฟฝ ๐ฌ๐ฌฬ ๐ฌ๐ฌฬ ๐๐ฬ ๏ฟฝ๏ฟฝ ๐ฎ๐ฎ๏ฟฝ d ๐ค๐ค ๐๐๐บ๐บ ๐ข๐ข โ ๏ฟฝ ๏ฟฝ๐๐ T ( ๐ฏ๐ฏ ) ๐๐ ฮผฮบ โ๐๐ T ( ๐ฏ๐ฏ ) ๐๐ ๏ฟฝ ๐๐๐บ๐บ ๐๐ ๐ง๐ง๏ฟฝ ๐ง๐ง๏ฟฝ ๐๐ฬ๐๐๐ค๐ค โ ๏ฟฝ ๏ฟฝ ๐๐ T ( ๐ฏ๐ฏ ) ๐๐ + ๐๐ T ( ๐ฏ๐ฏ ) ๐๐ ๏ฟฝ ๐ง๐ง๐๐๏ฟฝ d ๐ค๐ค ๐๐๐บ๐บ ๐๐ + ๏ฟฝ ๐๐ โ ๐๐ ๐๐๐บ๐บ ๐ข๐ข ๐ฏ๐ฏ T ๐ฎ๐ฎ๏ฟฝ d ๐ค๐ค + ๏ฟฝ ๐๐ ๐๐๐บ๐บ ๐๐ โ ๐๐ ๐๐๏ฟฝ T ( ๐ฏ๐ฏ ) ๐ง๐ง๏ฟฝ ๐๐ฬ d ๐ค๐ค , (31) โ ๏ฟฝ ๐๐ T ( ๐๐ ) ๐๐๏ฟฝ๐๐ ( ๐๐ ) ๐บ๐บ d ๐บ๐บ โ ๏ฟฝ ๐๐ T ( ๐๐ ) ๏ฟฝ๐๐ T ๐๐ ( ๐ฎ๐ฎ ) + ๐๐ ๐๐ ( ๐ฎ๐ฎ ) ๏ฟฝ ๐บ๐บ d ๐บ๐บ + ๏ฟฝ { ๐๐ T ( ๐๐ )} ๐๐ T ๐ง๐ง๏ฟฝ โฆ๐ฎ๐ฎโง d ๐ค๐ค ๐ค๐ค โ โช๐๐๐บ๐บ ๐ข๐ข โ ๏ฟฝ ๏ฟฝ ๏ฟฝ๐๐ , ( ๐๐ ) ๏ฟฝ๐๐ ๐ง๐ง๏ฟฝ โฆ๐ฎ๐ฎโง + ๏ฟฝ๐๐ , ( ๐๐ ) ๏ฟฝ๐๐ ๐ง๐ง๏ฟฝ โฆ๐ฎ๐ฎโง๏ฟฝ d ๐ค๐ค ๐ค๐ค โ โช๐๐๐บ๐บ ๐ข๐ข โ ๏ฟฝ ๏ฟฝ โฆ๐๐โง๐ง๐ง ๐๐T ๐๐๏ฟฝ { ๐๐ ( ๐๐ )}+{ ๐๐ T ( ๐๐ )} ๐๐๏ฟฝ๐ง๐ง ๐๐ โฆ๐๐โง๏ฟฝ d ๐ค๐ค ๐ค๐ค โ โช๐๐๐บ๐บ ๐๐ โ ๏ฟฝ โฆ๐๐โง๐ง๐ง ๐๐T ๏ฟฝ ๐๐ T { ๐๐ ( ๐ฎ๐ฎ )}+ ๐๐ { ๐๐ ( ๐ฎ๐ฎ )} ๏ฟฝ d ๐ค๐ค ๐ค๐ค โ โช๐๐๐บ๐บ ๐๐ + ๏ฟฝ { ๐๐ T ( ๐๐ )} ๐๐ ๐ง๐ง๏ฟฝ โฆ๐๐๏ฟฝ ( ๐ฎ๐ฎ ) โง ๐ค๐ค โ d ๐ค๐ค โ ๏ฟฝ ๏ฟฝ ๐๐ , ( ๐๐ ) ๐๐ ๐ง๐ง๏ฟฝ ๐ฌ๐ฌฬ ๐ฌ๐ฌฬ ๐๐ฬ + ๐๐ , ( ๐๐ ) ๐๐ ๐ง๐ง๏ฟฝ ๐ฌ๐ฌฬ ๐ฌ๐ฌฬ ๐๐ฬ ๏ฟฝ ๐ฎ๐ฎ d ๐ค๐ค ๐๐๐บ๐บ ๐ข๐ข + ๏ฟฝ ๐๐ T ( ๐๐ ) ๐๐ ๐ง๐ง๏ฟฝ ๐ง๐ง๏ฟฝ ๐ง๐ง๏ฟฝ ๐๐๏ฟฝ ( ๐ฎ๐ฎ ) ๐๐๐บ๐บ ๐๐ d ๐ค๐ค + ๏ฟฝ ๐๐ โ ๐๐ ๐๐๐๐ d ๐ค๐ค ๐๐๐บ๐บ ๐๐ + ๏ฟฝ ๐๐ โ ๐๐ โฆ๐๐โงโฆ๐๐โง d ๐ค๐ค ๐ค๐ค โ = โ ๏ฟฝ ๐๐๐๐ ๐บ๐บ d ๐บ๐บ โ ๏ฟฝ ๐๐๐๐๏ฟฝ ๐๐๐บ๐บ ๐๐ d ๐ค๐ค + ๏ฟฝ ๐๐ T ( ๐๐ ) ๐๐ T ๐ง๐ง๏ฟฝ ๐ฎ๐ฎ๏ฟฝ d ๐ค๐ค ๐๐๐บ๐บ ๐ข๐ข โ ๏ฟฝ ๏ฟฝ ๐๐ , ( ๐๐ ) ๐๐ ๏ฟฝ๐ง๐ง๏ฟฝ + ๐ง๐ง๏ฟฝ ๐ฌ๐ฌฬ ๐ฌ๐ฌฬ ๐๐ฬ ๏ฟฝ + ๐๐ , ( ๐๐ ) ๐๐ ๏ฟฝ๐ง๐ง๏ฟฝ + ๐ง๐ง๏ฟฝ ๐ฌ๐ฌฬ ๐ฌ๐ฌฬ ๐๐ฬ ๏ฟฝ๏ฟฝ ๐ฎ๐ฎ๏ฟฝ d ๐ค๐ค ๐๐๐บ๐บ ๐ข๐ข + ๏ฟฝ ๐๐ T ( ๐๐ ) ๐๐ ๐ง๐ง๏ฟฝ ๐ง๐ง๏ฟฝ ๐๐ฬ ๐๐๐บ๐บ ๐๐ d ๐ค๐ค โ ๏ฟฝ ๐๐ T ( ๐๐ ) ๐๐๏ฟฝ๐ง๐ง๐๐๏ฟฝ d ๐ค๐ค ๐๐๐บ๐บ ๐๐ + ๏ฟฝ ๐๐ โ ๐๐ ๐๐๐๐๏ฟฝ d ๐ค๐ค ๐๐๐บ๐บ ๐๐ . (32) where the normal and constant matrices ๐ง๐ง๏ฟฝ ๐๐ , ๐๐ฬ ๐๐ are shown in Appendix B. โ ๐๐ is a boundary-dependent parameter with the unit of length. Here it is defined as the distance between the points in subdomains sharing the boundary for ๐๐ โ ๐ค๐ค โ , and the smallest distance between the internal point and the external boundary for ๐๐ โ ๐๐ฮฉ . The penalty parameters are independent of the boundary size, in which ( ๐๐ , ๐๐ , ๐๐ , ๐๐ ) have the same unit of the Young's modulus, and ( ๐๐ , ๐๐ ) have the same unit of the permittivity of the material. As has been stated, the method is only stable when the penalty parameter ๐๐ is large enough. However, unlike the DG methods, in which the shape functions are completely independent, a small penalty parameter ๐๐ is enough to stabilize the present FPM. On the other hand, the accuracy may decrease if ๐๐ is too large. There is no upper limit for ๐๐ which corresponds to the essential boundary conditions. At last, the formula of the FPM can be converted into the following form: ๏ฟฝ๐๐ ๐ข๐ข๐ข๐ข ๐๐ ๐ข๐ข๐๐ ๐๐ ๐ข๐ข๐๐T ๐๐ ๐๐๐๐ ๏ฟฝ ๏ฟฝ๐ฎ๐ฎ๏ฟฝ๐๐๏ฟฝ ๏ฟฝ = ๏ฟฝ๐๐ ๐ข๐ข ๐๐ ๐๐ ๏ฟฝ , or ๐๐๐ฑ๐ฑ๏ฟฝ = ๐๐ , (33) where ( ๐ฎ๐ฎ๏ฟฝ , ๐๐๏ฟฝ ) are the global nodal displacement and electric potential vectors. ๐๐ is the global stiffness matrix, assembled by a series of point stiffness matrix ( ๐๐ ๐ธ๐ธ ) and boundary stiffness matrices ( ๐๐ โ , ๐๐ ๐ข๐ข , ๐๐ ๐๐ , ๐๐ ๐๐ , ๐๐ ๐ ๐ , ๐๐ ๐๐ and ๐๐ ๐๐ corresponding to ๐ค๐ค โ , ๐๐๐บ๐บ ๐ข๐ข , ๐๐๐บ๐บ ๐๐ , ๐๐๐บ๐บ ๐๐ , ๐๐๐บ๐บ ๐ ๐ , ๐๐๐บ๐บ ๐๐ and ๐๐๐บ๐บ ๐๐ respectively). ๐๐ is the global load vector. When establishing these stiffness matrices and load vectors, we use Gaussian quadrature rule. integration points are used in each subdomain, while only one integration point is used on each internal or external boundary. The numerical implementation is omitted here. The final stiffness matrix is sparse and symmetric, which is good for modeling complex problems with a large number of DOFs.
5. NUMERICAL EXAMPLES 5.1 Hollow cylinder
The first example is an infinite length flexoelectric tube. As shown in Fig. 3, this is a plane strain problem with axisymmetric boundary conditions. The material properties are: Young's modulus ๐ธ๐ธ =139GPa , Poisson's ratio ๐๐ = 0.3 , internal material length ๐๐ =2 ๐๐๐๐ , flexoelectric coefficients ๐๐ฬ = ๐๐ฬ = ๐๐ฬ = 1 ร10 โ6 C/m and permittivity of the dielectric ๐ ๐ ฬ = ๐ ๐ ฬ =1 ร 10 โ9 F/m (see Appendix A). The material is non- piezoelectric. The geometric parameters are: ๐๐ ๐๐ = 10 ๐๐๐๐ , ๐๐ ๐๐ =20 ๐๐๐๐ . And the boundary conditions are given as: ๐ข๐ข ๐๐ =0.045 ๐๐๐๐ , ๐ข๐ข ๐๐ = 0.05 ๐๐๐๐ , ๐๐ ๐๐ = 0 ๐๐ , ๐๐ ๐๐ = 1 ๐๐ . In the FPM, according to the symmetry, only a quarter of the entire domain is considered. Symmetric boundary conditions are applied. 1260 uniform points are employed. The computational parameters are given as: ๐๐ = โ , ๐๐ = 1 ร 10 ๐ธ๐ธ , ๐๐ =1 ร 10 ๐ ๐ ฬ , ๐๐ = 2.0 ๐ธ๐ธ , ๐๐ = 100 ๐ธ๐ธ , ๐๐ = 0 . 8 FIGURE 3:
AN INFINITE LENGTH TUBE WITH AN AXISYMMETRIC CROSS SECTION. (A) (B)
FIGURE 4:
THE COMPUTED SOLUTION. (A) DISTRIBUTION OF RADIAL DISPLACEMENT. (B) DISTRIBUTION OF ELECTRIC POTENTIAL.
The same example is also considered in [2] and [4] using the mixed FEM. The computed solution by the FPM is presented in Fig. 4, comparing with the analytical solution given by Mao and Purohit [17]. As can be seen, the numerical solution shows excellent agreement with the analytical result. We define the relative errors of displacement and electric potential as: ๐๐ ๐ข๐ข = โ๐ฎ๐ฎ โ โ ๐ฎ๐ฎโ ๐ฟ๐ฟ โ๐ฎ๐ฎโ ๐ฟ๐ฟ , ๐๐ ๐๐ = โ๐๐ โ โ ๐๐โ ๐ฟ๐ฟ โ๐๐โ ๐ฟ๐ฟ , (34) where โ๐ฎ๐ฎโ ๐ฟ๐ฟ = ๏ฟฝ ๐ฎ๐ฎ T ๐ฎ๐ฎ d ๐บ๐บ ๐บ๐บ , โ๐๐โ ๐ฟ๐ฟ = ๏ฟฝ ๐๐ d ๐บ๐บ ๐บ๐บ In this example, the relative errors are ๐๐ ๐ข๐ข = 3.2 ร 10 โ5 and ๐๐ ๐๐ = 8.4 ร 10 โ4 . Apparently, the present FPM presents significantly high accuracy in this benchmark example. FIGURE 5:
2D BLOCK SUBJECT TO A POINT LOAD. (A) (B)
FIGURE 6:
THE COMPUTED SOLUTION. (A) DISTRIBUTION OF ELECTRIC POTENTIAL ๐๐ ( ๐๐ ). (B) DISTRIBUTION OF ELECTRIC FIELD ๐ธ๐ธ ( ๐๐ ๐๐๐๐โ ). Second, we consider a 2D plane strain model of a block subject to a concentrated load. As shown in Fig. 5, a significant strain gradient and electric field can be expected. The example is initiated from the flexoelectric effect observed in some atomic force microscope experiments [2].
The material properties are the same as the first example. The geometric parameters are: ๐๐ = 20 ๐๐๐๐ , ๐๐ = 10 ๐๐๐๐ . The concentrated force ๐น๐น = 100 ๐๐๐๐ , applied uniformly in a ๐๐๐๐ width area to avoid singularity. 3200 points are utilized in the 9 FPM. Both Voronoi Diagram partition and quadrilateral partition (generated by the preprocessing module of ABAQUS) are employed. The result shows great consistency between the different partitions. For simplicity, only the solution of the quadrilateral partition is presented in the paper. The computational parameters in the FPM are: ๐๐ = โ , ๐๐ = 1 ร 10 ๐ธ๐ธ , ๐๐ = 1 ร 10 ๐ ๐ ฬ , ๐๐ = 1.0 ๐ธ๐ธ , ๐๐ =50 ๐ธ๐ธ , ๐๐ = 0 . Figure 6 exhibits the numerical solution of the electric potential ๐๐ and the vertical electric field ๐ธ๐ธ . The result shows great consistency with [2] in which the example is studied by the mixed FEM. The third example is a truncated pyramid subject to a uniform force (see Fig. 7). In this example, the material is anisotropic. And the piezoelectric effect is taken into consideration. The material properties are: ๐ธ๐ธ = 100GPa , ๐๐ =0.37 , ๐ ๐ ฬ = 11 ร 10 โ9 F/m , ๐ ๐ ฬ = 12.48 ร 10 โ9 F/m , ๐๐ฬ = 1 ร 10 โ6 C/m , ๐๐ฬ = ๐๐ฬ = 0 and piezoelectric coefficients ๐๐ = โ โ , ๐๐ = ๐๐ = 0 . The strain gradient effect is neglected ( ๐๐ = 0 ). The geometric parameters are: ๐๐ = 750 ๐๐๐๐ , ๐๐ = 2250 ๐๐๐๐ , ๐๐ = 750 ๐๐๐๐ . The force ๐น๐น = 450 ๐๐๐๐ is distributed uniformly on the top and the bottom surface. The electric potential on the top surface is fixed to zero while on the bottom surface, it is subject to a constant but priori unknown value ๐๐ . FIGURE 7:
SIMPLY SUPPORTED TRUNCATED PYRAMID.
With 306 uniform points and quadrilateral partition used, the numerical solutions of the electric potential ๐๐ and mechanical strain ๐๐ are presented in Fig. 8. The computational parameters ๐๐ = โ , ๐๐ = 1 ร 10 ๐ธ๐ธ , ๐๐ = 1 ร 10 ๐ ๐ ฬ , ๐๐ =1.0 ๐ธ๐ธ , ๐๐ = 0 , ๐๐ = ๐ ๐ ฬ . The deformation of the truncated pyramid shows a bending component as expected. The mechanical strain gradient is significant, as well as the nonhomogeneous distribution of the electric potential. The same example is also studied in [2, 10] based on LME meshfree approximation, in which the numerical result agrees well with the current result achieved by the FPM. (A) (B) FIGURE 8:
THE COMPUTED SOLUTION. (A) DISTRIBUTION OF ELECTRIC POTENTIAL ๐๐ ( ๐๐ ). (B) DISTRIBUTION OF STRAIN ๐๐ .
6. CONCLUSION
A new computational approach is developed for analyzing piezoelectric and flexoelectric effects in 2D solids. The purely meshless Fragile Points Method (FPM) based on Galerkin weak-form formulation is superior to the previous meshless methods (e.g., EFG, MLPG, etc.) as it has a much simpler numerical integration scheme based on a polynomial shape function, and results in symmetric and sparse matrices. Numerical examples are carried out for validation. In conclusion, the present FPM is accurate, efficient, and has a great potential in analyzing flexoelectric problems in complex geometries.
REFERENCES [1] Wang, Bo, Gu, Yijia, Zhang, Shujun, and Chen, Long-Qing. โFlexoelectricity in solids: Progress, challenges, and perspectives.โ
Progress in Materials Science
Vol. 106 (2019): pp. 100570. [2] Zhuang, Xiaoying, Nguyen, Binh Huy, Nanthakumar, Subbiah Srivilliputtur, Tran, Thai Quoc, Alajlan, Naif, and Rabczuk, Timon. โComputational Modeling of FlexoelectricityโA Review.โ
Energies
Vol. 13 No. 6 (2020): pp. 1326. [3] Mao, Sheng, Purohit, Prashant K., and Aravas, Nikolaos. โMixed finite-element formulations in piezoelectricity and flexoelectricity.โ
Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences
Vol. 472 (2016): pp. 20150879. [4] Deng, Feng, Deng, Qian, Yu, Wenshan, and Shen, Shengping. โMixed finite elements for flexoelectric solids.โ
Journal of Applied Mechanics
Vol. 84 No. 8 (2017): pp. 081004. [5] Maranganti, R., Sharma, N. D., and Sharma, P.. โElectromechanical coupling in nonpiezoelectric materials due to nanoscale nonlocal size effects: Greenโs function solutions and embedded inclusions.โ
Physical Review B
Vol. 74 No. 1 (2006): pp. 014110. [6] Hu, ShuLing, and Shen, ShengPing. โVariational principles and governing equations in nano-dielectrics with the 10 flexoelectric effect.โ
Science China Physics, Mechanics and Astronomy
Vol. 53 No. 8 (2010): pp. 1497-1504. [7] Deng, Qian, Liu, Liping, and Sharma, Pradeep. โFlexoelectricity in soft materials and biological membranes.โ
Journal of the Mechanics and Physics of Solids
Vol. 62 (2014): pp. 209-227. [8] Yvonnet, Julien, and Liu, L. P.. โA numerical framework for modeling flexoelectricity and Maxwell stress in soft dielectrics at finite strains.โ
Computer Methods in Applied Mechanics and Engineering
Vol. 313 (2017): pp. 450-482. [9] Sladek, Jan, Sladek, Vladimir, Wรผnsche, Michael, and Zhang, Chuanzeng. โEffects of electric field and strain gradients on cracks in piezoelectric solids.โ
European Journal of Mechanics-A/Solids
Vol. 71 (2018): pp.187-198. [10] Abdollahi, Amir, Peco, Christian, Millan, Daniel, Arroyo, Marino, and Arias, Irene. โComputational evaluation of the flexoelectric effect in dielectric solids.โ
Journal of Applied Physics
Vol. 116 No. 9 (2014): pp. 093502. [11] He, Bo, Javvaji, Brahmanandam, and Zhuang, Xiaoying. โCharacterizing Flexoelectricity in Composite Material Using the Element-Free Galerkin Method.โ
Energies
Vol. 12 No. 2 (2019): pp. 271. [12] Tang, Z., Shen, S., and Atluri, S. N.. โAnalysis of materials with strain-gradient effects: A meshless local Petrov-Galerkin (MLPG) approach, with nodal displacements only.โ
Computer Modeling in Engineering and Sciences
Vol. 4 No. 1 (2003): pp. 177-196. [13] Hillman, Michael, and Chen,
JiunโShyan. โAn accelerated, convergent, and stable nodal integration in Galerkin meshfree methods for linear and nonlinear mechanics.โ
International Journal for Numerical Methods in Engineering
Vol. 107 No. 7 (2016): pp. 603-630. [14] Duan, Qinglin, Li, Xikui, Zhang, Hongwu, and Belytschko, Ted. โ
Secondโorder accurate derivatives and integration schemes for meshfree methods.โ
International Journal for Numerical Methods in Engineering
Vol. 92 No. 4 (2012): pp. 399-424. [16] Shu, C., Ding, H., and Yeo, K. S.. โLocal radial basis function-based differential quadrature method and its application to solve two-dimensional incompressible NavierโStokes equations.โ
Computer methods in applied mechanics and engineering
Vol. 192 No. 7-8 (2003): pp. 941-954. [17] Mozolevski, Igor, Sรผli, Endre, and Bรถsing, Paulo R. โhp-version a priori error analysis of interior penalty discontinuous Galerkin finite element approximations to the biharmonic equation.โ
Journal of Scientific Computing
Vol. 30 No. 3 (2007): pp. 465-491. [16] Mao, Sheng, and Purohit, Prashant K.. โInsights into flexoelectric solids from strain-gradient elasticity.โ
Journal of Applied Mechanics
Vol. 81 No. 8 (2014): pp. 081004.
Appendix A: Material Properties
The material property matrices are given as: D ๏ฟฝ ฯฮต = ๏ฟฝ๐๐ + 2 ๐บ๐บ ๐๐ ๐๐ ๐๐ + 2 ๐บ๐บ
00 0
๐บ๐บ๏ฟฝ , D ๏ฟฝ ฮผฮบ = ๐๐ โโฃโขโขโขโขโขโขโขโขโก ๐๐ + 2 ๐บ๐บ ๐๐ ๐๐ + 2 ๐บ๐บ ๐๐ ๐๐ ๐๐ +3 ๐บ๐บ ๐บ๐บ ๐๐ ๐๐ +3 ๐บ๐บ ๐บ๐บ
00 0 0 ๐บ๐บ ๐บ๐บ
00 0 ๐บ๐บ ๐บ๐บ โฆโฅโฅโฅโฅโฅโฅโฅโฅโค , ๐๐๏ฟฝ = ๏ฟฝ๐ ๐ ฬ ๐ ๐ ฬ ๏ฟฝ , ๐๐ = ๏ฟฝ ๐๐ ๐๐ ๐๐ ๏ฟฝ , (35) ๐๐ = ๏ฟฝ๐๐ฬ ๐๐๏ฟฝ +๐๐๏ฟฝ ๐๐ฬ ๐๐ฬ
11 ๐๐๏ฟฝ +๐๐๏ฟฝ ๐๐ฬ ๏ฟฝ T , where ( ๐๐ , ๐บ๐บ ) are Lamรฉ parameters, ๐๐ is the internal material length, ๐๐๏ฟฝ , ๐๐ , and ๐๐ are related to the permittivity of the dielectric, the piezoelectric tensor and the flexoelectric tensor respectively. Appendix B: Matrices in Numerical Implementation ๐๐ฬ = ๏ฟฝ ๏ฟฝ , ๐๐ฬ = ๏ฟฝ ๏ฟฝ , ๐๐ฬ = ๏ฟฝ โ โ ๏ฟฝ , ๐๐ฬ = ๏ฟฝ โ โ ๏ฟฝ , ๐ง๐ง๏ฟฝ = ๏ฟฝ๐๐ ๐๐ ๐๐ ๐๐ ๏ฟฝ , (36) ๐ง๐ง๏ฟฝ = ๏ฟฝ๐๐ ๐๐ ๐๐ ๐๐ ๏ฟฝ , ๐ง๐ง๏ฟฝ = ๏ฟฝ ๐๐ ๐๐ ๐๐ ๐๐ ๏ฟฝ , ๐ง๐ง๏ฟฝ = ๏ฟฝ๐๐ ๐๐ ๐๐ ๐๐ ๐๐ ๐๐
00 0 0 ๐๐ ๐๐ ๏ฟฝ , ๐ง๐ง๏ฟฝ = ๏ฟฝ๐๐ ๐๐ ๐๐ ๐๐ ๏ฟฝ , ๐ฌ๐ฌฬ = ๏ฟฝ๐๐ โ๐๐ โ๐๐ ๐๐ ๏ฟฝ ..