A New Meshless Fragile Points Method (FPM) With Minimum Unknowns at Each Point, For Flexoelectric Analysis Under Two Theories with Crack Propagation. Part I: Theory and Implementation
aa r X i v : . [ m a t h . NA ] S e p A new meshless fragile points method (FPM) with minimumunknowns at each point, for flexoelectric analysis under two theorieswith crack propagation. Part I: Theory and implementation
Yue Guan a, ∗ , Leiting Dong b , Satya N. Atluri a a Department of Mechanical Engineering, Texas Tech University, Lubbock, TX 79415, United States b School of Aeronautic Science and Engineering, Beihang University, Beijing 100191, China
Abstract “Flexoelectricity” refers to a phenomenon which involves a coupling of the mechanical strain gra-dient and electric polarization. In this study, a meshless Fragile Points Method (FPM), is presentedfor analyzing flexoelectric effects in dielectric solids. Local, simple, polynomial and discontinu-ous trial and test functions are generated with the help of a local meshless Differential Quadratureapproximation of derivatives. Both primal and mixed FPM are developed, based on two alternateflexoelectric theories, with or without the electric gradient effect and Maxwell stress. The first the-ory is fully nonlinear and is recommended at the nano-scale, while the second theory is linear andis sufficient at the micro-scale. In the present primal as well as mixed FPM, only the displacementsand electric potential are retained as explicit unknown variables at each internal Fragile Point inthe final algebraic equations. Thus the number of unknowns in the final system of algebraic equa-tions is kept to be absolutely minimal. An algorithm for simulating crack initiation and propagationusing the present FPM is presented, with classic stress-based criterion as well as a Bonding-Energy-Rate(BER)-based criterion for crack development. The present primal and mixed FPM approachesrepresent clear advantages as compared to the current methods for computational flexoelectric analy-ses, using primal as well as mixed Finite Element Methods, Element Free Galerkin (EFG) Methods,Meshless Local Petrov Galerkin (MLPG) Methods, and Isogeometric Analysis (IGA) Methods, be-cause of the following new features: they are simpler Galerkin meshless methods using polynomialtrial and test functions; minimal DoFs per Point make it very user-friendly; arbitrary polygonal sub-domains make it flexible for modeling complex geometries; the numerical integration of the primal ∗ Corresponding author.
Email address: [email protected] (Yue Guan)
Preprint submitted to November 13, 2020 s well as mixed FPM weak forms is trivially simple; and FPM can be easily employed in crackdevelopment simulations without remeshing or trial function enhancement. In this first part of thetwo-part series, we focus on the theoretical formulation and implementation of the proposed primalas well as mixed FPM. Numerical results and validation are then presented in Part II of the presentpaper.
Keywords:
Flexoelectricity, Strain gradient effect, Fragile Points Method (FPM), Crack propagation
1. Introduction “Flexoelectricity” describes an electromechanical coupling effect between the electric polariza-tion and mechanical strain gradient [1–3]. As a result of the recent trend of miniaturization ofelectromechanical systems, which necessitates the consideration of the significant effects of straingradients, the flexoelectric effect has gained much attention in recent years. Resulting from thesize effect, the flexoelectric response can become significant and even dominant in micro/nano-electromechanical systems [4–6]. Therefore, there has been an increased demand for reliable andaccurate theories and numerical methods for flexoelectric analyses. In 2006, Maranganti et al. [3]proposed the first continuum theory for flexoelectricity. After that, a number of advanced contin-uum theories have been developed, considering the Maxwell stress and surface effect [7, 8] andincorporating kinetic energy [9], etc.These continuum theories of flexoelectricity have been applied in analyzing mechanical andelectrical responses of various systems. For example, Mao and Purohit [10] derived the governingequations for a flexoelectric solid under small deformation and presented the analytical solutionsfor a 2D axisymmetric boundary value problem. The electric field gradient effect, as well as theconverse flexoelectricity (a linear coupling effect between the electric field gradient and mechanicalstrain) are omitted in the governing equations. In the present paper, we denote these equations in [10]as a “reduced” theory. On the contrary, the generalized form considering both the electric gradienteffect and the electroelastic stress is signified as a “full” theory [2, 3]. Whereas the reduced theoryis linear and much easier for numerical implementation, the influence of the electric gradient effectand the electroelastic stress can be significant at the nano-scale and thus the reduced theory maylead to inaccurate results [8]. In this paper, both the nonlinear and linear theories are numerically2mplemented. The appropriate theory should be chosen depending on the length scale of the problemunder consideration.In practice, the governing equations and boundary value problems describing the flexoelectriceffect can be formulated from two different approaches. First, as presented by Maranganti et al. [3],the internal energy density can be given as a function of strain, strain gradient, electric polarizationand polarization gradient. Since the electric polarization ( P ) is used as an independent variable,this is denoted as a ‘P-formulation’. Alternatively, we can employ the electric Gibbs free energyto formulate the boundary value problem. Thus instead of the electric polarization and polariza-tion gradient, the independent variables become electric field ( E ) and electric field gradient, andthe resulting equations are known as an âĂŸE-formulationâĂŹ. Both the formulations have beenconsidered for numerical implementations. For instance, Yvonnet and Liu [11] and Mao et al. [12]have proposed primal and mixed finite element (FEM) frameworks based on the P- formulation.Whereas the E-formulation can also be incorporated with multiple numerical discretization meth-ods, including primal and mixed FEM [13, 14], Element-Free-Galerkin (EFG) method [15] andMeshless Local Petrov-Galerkin (MLPG) method [16], etc. The two formulations are inherentlyequivalent and can be unified by a linear relation between the electric polarization and electric dis-placement (the conjugate component of the electric field). Whereas the P-formulation employs theelectric polarization induced by strain gradient intuitively, when the full theory is taken into con-sideration, the E-formulation is more favorable due to mathematical simplicity [2].. In the currentwork, we also concentrate on the E-formulation.As have been stated, a number of numerical methods have been developed in analyzing flex-oelectric and / or strain gradient effects in dielectric solids. These methods can be divided intothree groups: the traditional Finite Element Methods (FEM), meshless methods, and isogeometricanalysis (IGA) methods. Generally, the flexoelectric behavior is governed by a fourth-order partialdifferential equation (PDE). As a result, the main difficulty in modelling flexoelectricity for classicelement-based methods lies in the C continuity requirement. On one hand, the C requirementcan be satisfied by using Argyris triangular element [11] or subparametric quadrilateral element[17, 18]. However, even for two dimensional problems, the C continuous elements result in largenumbers of degrees of freedom (DoFs); for both of the approaches [11, 17], 9 DoFs are required ateach node, including 6 mechanical quantities and 3 electrical quantities, thus making the formula-3ion rather unfriendly for users. Besides, generating C continuous elements in 3D can be extremelyintractable and has not been pursued in any prior literature. Alternatively, mixed FEM can be devel-oped using displacement gradients and Lagrangian multipliers as additional independent variablesand requiring only C continuity of all the variables. These mixed FEM have been applied in an-alyzing strain gradient [19–21], flexoelectric [12, 22] and converse flexoelectric effects [23], andhave also been extended to 3D analysis [14] and topology optimization [24]. Yet they have evenmore DoFs in each element (e.g., 87 DoFs in a 9-node element [12], 36 DoFs in a 5-node element[23], or 37 DoFs in a 7-node element [22]), making it prohibitively complex for practical use.Another category of methods, known as “meshless methods”, is partly or completely free ofmesh discretization and have also shown their capability in analyzing flexoelectric or strain-gradienteffects. Abdollahi et al. [25] proposed a meshfree method with a C ∞ continuous basis function us-ing Local Maximum Entropy (LME) approximants. After that, the Element-Free Galerkin (EFG)Method based on the Moving Least Squares (MLS) approximation which can generate C con-tinuos functions through an appropriate choice of the weight functions, was applied to compositebeam analysis with flexoelectricity [26–28], as well as 2D strain gradient [29] and flexoelectric anal-ysis [15]. The Meshless Local Petrov-Galerkin (MLPG) Method is also established based on theMLS approximation for trial functions, whereas multiple test functions (e.g., weight function, lo-cal fundamental solution, Heaviside function, etc.) can be applied and thus lead to an asymmetricformulation. The MLPG method is also easily applicable in analyzing strain gradient [30] and flex-oelectric [16, 31] problems. The main advantage of the MLS approximation used in the EFG andMLPG methods is with an appropriate weight function, we can generate smooth trial functions with C or higher order continuity. However, the trial functions are very complicated rational polynomi-als, and its higher order derivatives can be rational, and even more complicated. In order to avoidthat problem, a mixed MLPG method, or Meshless Finite Volume Method (FVM) is proposed byAtluri et al. [32], with MLS trial functions and Heaviside test functions. Several other kinds ofmixed MLPG methods were also developed for solving fourth order ODEs or PDEs arising in for-mulations involving strain gradient effects [33, 34]. In these mixed MLPG methods, even thoughthe mechanical strain and / or strain gradients are used as independent variables, these additionalvariables are eliminated at the node level and only the displacement DoFs are retained at the nodesin the final formulations. Nevertheless, even though higher order derivatives are avoided, the trial4unction based on MLS approximation is still complicated in the mixed MLPG method, making itdifficult for high-precision numerical integration. Similar argument can be made when CompactlySupported Radial Basis Function (CSRBF) is used as the trial function when solving nonlinearflexoelectric problems [35].All the meshless methods mentioned above have certain drawbacks: First, the shape functionsbased on MLS, LME lack the Kronecker delta property. Therefore, additional treatments are re-quired in imposing the essential boundary conditions, e.g., a modified collocation method based oninterior penalty functions [36] or Lagrange multipliers [28]. The imposition of higher order essentialboundary conditions in flexoelectric analysis has not been mentioned in any meshless approaches.Second, as a result of the complicated trial functions, the numerical integration of the weak formin either EFG or MLPG methods can be tedious. Actually, this difficulty in domain integration is achallenge for most Galerkin meshfree methods [37]. The simplest choice is direct nodal integration.Though efficient, and no background mesh required, the direct nodal integration may have stabilityissues and can be less accurate or non-convergent. High order quadrature, on the other hand, canachieve stability and better convergence. But it is computationally expensive and prohibitive to beused in practice. A number of modified nodal integration methods are then proposed, in order toensure the accuracy and stability, which in turn sacrifice the efficiency [37]. Some other newer-typequadrature, e.g., a three-point integration scheme using background triangle elements [38], are alsoadopted. Nevertheless, these techniques are all applied to remedy the complicated integration in theweak form for the MLS trial functions.Alternatively, the difficulty of numerical integration can be simply solved by employing simpleshape functions like polynomials. The Fragile Points Method (FPM), in contrast to all the previousmeshless methods, is based on simple, polynomial, piecewise-continuous trial and test functions[39, 40]. Consequently, the classic Gaussian quadrature can be employed, and numerical integra-tions are trivially simple. Only one integration point in each subdomain is sufficient most of thetime in the FPM, whereas in methods like EFG and MLPG, even a large number of Gaussian inte-gration points cannot guarantee a reliable numerical integration. Furthermore, the shape functionsin the FPM pass through the nodal data, thus the Dirichlet boundary conditions can be enforceddirectly. On the other hand, when compared with the element-based methods, the trial functionsof the FPM are generated using scattered points with the support of each point and always involve5he same Cartesian strains in each subdomain. Thus, it has benefits in avoiding the mesh distor-tion and locking problems associated with element-based methods. Besides, unlike the FEM andother element-based methods which require triangular or quadrilateral element meshing, the FPMcan be implemented with random point distributions and arbitrary subdomain shapes, and thus candeal with problems with complicated overall geometries. The FPM has already shown great ac-curacy and efficiency in solving 2D heat conduction [40, 41] and elasticity problems [39, 42]. Inthis study, we formulate and apply the FPM for analyzing flexoelectric problems using both full andreduced theories. Furthermore, we can also develop mixed formulations for the FPM in which themechanical strain and strain gradients are used as independent variables. As in the mixed MLPGmethod, the additional higher order variables can all be eliminated at each Point level. Therefore,the final mixed FPM formulations simply have nodal displacement and electric potentials, i.e., only3 DoFs are retained at each node for a 2D problem in the final analysis. Comparisons of all therepresentative primal and mixed numerical methods are shown in Tables 1 and 2 respectively.Otherwise, the flexoelectric effect can also be analyzed by the IGA approach [6, 43], in whichNon-Uniform Rational B-Spline (NURBS) basis functions with C or higher order continuity arechosen to approximate both the geometry and field variables (displacement and electric potentialfields). Yet the NURBS basis function also lacks Kronecker-delta property and thus, as in EFG,MLPG, etc., the essential boundary conditions need to be imposed by Lagrange multipliers [44].The IGA approach also requires a large number of DoFs in each patch and thus making it ratherunfriendly for users.All the above studies just concentrate on continuous domains. Yet a significant strain gradientcan be expected in the near-tip fields for a crack in elastic materials. Thus the strain gradient andflexoelectric effects should be taken into consideration for fracture analysis at micro and nano scales.Sladek et al. has analyzed the flexoelectric effect for stationary cracks in dielectric solids with theFEM [18] and MLPG method [45]. Abdollahi and Arias has studied the crack propagation for afour-point bending test in piezoelectric material [46], as well as the influence of flexoelectricity ona stationary crack [47] based on a meshless formulation. However, few studies have been reportedon the crack propagation simulations considering the flexoelectric effect. This may stem from thedifficulties in simulating crack propagation while using continuous trial functions. Approaches tosimulate crack and rupture initiation and propagations, e.g. Generalized FEM (GFEM) [48], Ex-6ended FEM (XFEM) [49, 50], or Zencrack [51], usually involve mesh refinement or trial functionenrichment. Nevertheless, in the present FPM, as a result of the discontinuous trial and test func-tions, the algorithm for simulating crack propagation is relatively simple and intuitive, and does notinvolve either mesh refinement or trial function enrichment.In this paper, we focus on analyzing the flexoelectric effect in dielectric solids using the FPM.Both full and reduced flexoelectric theories are considered and their corresponding equations areintroduced in section 2. The primal and mixed FPM formulations are developed in sections 3 and4 respectively. Section 5 presents the algorithm for simulating crack initiation and propagationwhile using the FPM. Numerical results and validation of the overall methodologies, as well as aparametric discussion, are given in Part II of this 2-part paper.
2. Flexoelectricity theories and boundary value problems
In the current work, we consider a homogeneous elastic dielectric 2D domain Ω with bound-ary ∂ Ω . The spatial coordinate is given as x = [ x , x ] T . Subjected to mechanical and electricalloadings, the responses of the material are described by the displacement vector field u ( x , x ) =[ u , u ] T and the electric potential field φ ( x , x ) .Phenomenologically, the flexoelectric effect describes an electric polarization generated by themechanical strain gradient: P i = µ ijkl κ jkl . (1)Einstein summation convention is used here. P = [ P , P ] T is the electrical polarization vector, κ jkl = u j,kl is the second gradient of displacement (i.e., the strain gradient), and µ ijkl are elementsof a fourth order flexoelectric tensor ( i, j, k, l = 1 , ).The present boundary value problem is formulated from an electric Gibbs free energy G , inwhich the mechanical strain ε , strain gradient κ , electric field E and electric field gradient V are7 able 1: A comparison of primal numerical methods. Method Primal FPM Primal FEM Primal EFG PrimalMLPG IGATrial func-tions Discontinuouspolynomial C · · · As aboveDelta func-tion property Yes Yes Yes/No Yes/No NoWeak form Galerkin(presentpaper) orPetrov-Galerkin(future) Galerkin Galerkin Petrov-Galerkin GalerkinWeak formintegration Simple(Gaussianquadrature) Simple(Gaussianquadrature) Very compli-cated Very compli-cated Very compli-catedDoF at eachnode 3 9 6 3 38 able 2: A comparison of mixed numerical methods.
Method Mixed FPM Mixed FEM Mixed MLPGTrial functions Discontinuous linearpolynomial C linear/quadraticpolynomial MLS, etc.Test functions As above As above MLS, weight func-tion, local funda-mental solution,Dirac delta function,Heaviside function,etc.Delta functionproperty Yes Yes Yes/NoWeak form Galerkin (presentpaper) or Petrov-Galerkin (future) Galerkin Petrov-GalerkinWeak form inte-gration Very simple (one-point quadrature) Simple (Gaussianquadrature) Very complicatedDoF at eachnode 3 3 - 13 39ndependent variables: ε = [ ε , ε , ε ] T , κ = [ κ , κ , κ , κ , κ , κ ] T , E = [ E , E ] T , V = [ V , V , V , V ] T . ε ij = 12 ( u i,j + u j,i ) ,κ ijk = u i,jk ,E i = − φ ,i ,V ij = E i,j . (2)The governing equations are given as: (cid:0) σ ij − µ ijk,k + σ ESij (cid:1) ,j + b i = 0 , (3) ( D i − Q ij,j ) ,i = q, (4)where b i is the body force per volume, q is the free charge per volume, σ ij , µ ijk , D i and Q ij are the(local) stress, higher-order stress, electric displacement, and higher-order electric displacement, i.e.,conjugate variables of ε ij , κ ijk , E i and V ij respectively. σ ESij is the generalized electrostatic stress,which can be written as: σ ESij = σ Mij + µ Mijm,m ,σ Mij = E i D j − E k D k δ ij ,µ Mijm,m = V ik Q kj − V kl Q kl δ ij . (5)where σ Mij and µ Mijm,m are the Maxwell stress and the higher-order electrostatic stress respectively.The corresponding Dirichlet boundary conditions are: u i = e u i , on ∂ Ω u , (6) φ = e φ, on ∂ Ω φ , (7) ∇ n u i = u i,l n l = e d i , on ∂ Ω d , (8) ∇ n φ = φ ,l n l = e P , on ∂ Ω P . (9)The Neumann boundary conditions are: (cid:0) σ ij − µ ijk,k + σ ESij (cid:1) n j + ∇ tk ( n l ) n m n j µ ijm − ∇ tj ( n m µ ijm ) = e Q i , on ∂ Ω Q , (10) ( D i − Q ij,j ) n i + ∇ tl ( n l ) n i n j Q ij − ∇ ti ( Q ij n j ) = − e ω, on ∂ Ω ω , (11) µ ijm n m n j = e R i , on ∂ Ω R , (12)10 ij n i n j = e Z, on ∂ Ω Z . (13)In the above equations, e u , e φ , e d , e P , e Q , e ω , e R and e Z are known functions. n = [ n , n ] T is theoutward unit normal vector to ∂ Ω . ∇ n = n · ∇ and ∇ t = ∇ − n ∇ n are the normal and tangentialsurface differential operators respectively. ∂ Ω u ∪ ∂ Ω Q = ∂ Ω φ ∪ ∂ Ω ω = ∂ Ω d ∪ ∂ Ω R = ∂ Ω P ∪ ∂ Ω Z = ∂ Ω , and ∂ Ω u ∩ ∂ Ω Q = ∂ Ω φ ∩ ∂ Ω ω = ∂ Ω d ∩ ∂ Ω R = ∂ Ω P ∩ ∂ Ω Z = ∅ .The constitutive relations are given in the matrix forms: σ = C σε ε − eE − bV , µ = C µκ κ − aE , D = ΛE + e T ε + a T κ , Q = ΦV + b T ε . (14)where σ = [ σ , σ , σ ] T , µ = [ µ , µ , µ , µ , µ , µ ] T , D = [ D , D ] T , Q = [ Q , Q , Q , Q ] T . C σε , C µκ , Λ , Φ , a , b and e are matrices of material properties. The generalized electrostatic stress σ ES can also be written in a matrix form: σ ES = (cid:2) σ ES , σ ES , σ ES , σ ES (cid:3) T = b DE + b QV , (15)where b D = D − D − D D D D , b Q = Q − Q Q − Q − Q − Q Q − Q − Q
12 12 Q Q Q Q Q . Now, some different formulations can be reduced from the above-mentioned generalized formu-lation in the absence of electrostatic stress σ ES and / or of the electric field gradient V . Here we11resent a reduced formulation [10] in which both σ ES and V are omitted. The governing equationsand boundary conditions in the reduced form are: ( σ ij − µ ijk,k ) ,j + b i = 0 , (16) D i,i = q, (17)Dirichlet boundary conditions: u i = e u i , on ∂ Ω u , (18) φ = e φ, on ∂ Ω φ , (19) ∇ n u i = u i,l n l = e d i , on ∂ Ω d , (20)Neumann boundary conditions: (cid:0) σ ij − µ ijk,k + σ ESij (cid:1) n j + ∇ tk ( n l ) n m n j µ ijm − ∇ tj ( n m µ ijm ) = e Q i , on ∂ Ω Q , (21) D i n i = − e ω, on ∂ Ω ω , (22) µ ijm n m n j = e R i , on ∂ Ω R . (23)And the corresponding constitutive relations are: σ = C σε ε − eE , µ = C µκ κ − aE , D = ΛE + e T ε + a T κ . (24)As the electrostatic stress σ ES is omitted, this reduced form is a linear formulation, and thusis convenient for numerical implementation. It remains accurate in micro or larger scale structureswith negligible electric field gradient. However, studies have shown that the significance of theelectric force increases dramatically with the minimization of the system size [8]. Therefore, fornano-structures or structures with significant electric field gradient, the full theory including theelectrostatic stress should be taken into consideration.Note that an alternative formulation (“P-formulation”) based on an internal energy e U whichused mechanical strain ε , strain gradient κ , electrical polarization P and polarization field P ,i asindependent variables are also commonly used in previous literatures [3, 11]. Which formulation12eing used is simply a matter of choice, since the two formulations are essentially equivalent. Andthe relation between the electric displacement, polarization and field are given as: D = ǫ E + P , (25)where ǫ is the permittivity of free space.
3. The Primal Fragile Points Method (primal FPM)
In the Fragile Points Method (FPM), first, a set of random Points are scattered in the domain.The entire domain can then be partitioned into several nonoverlapping subdomains. Within eachsubdomain, only one Fragile Point exists. Multiple partitioning schemes may be employed, e.g.,the Voronoi Diagram partition (see Fig. 1(a)), as well as quadrilateral and triangular partition. Thetraditional FEM meshing can also be employed. The FEM elements can be converted into FPMsubdomains, while the internal Fragile Point is defined as the centroid of each FEM geometricalelement. Thus, the preprocessing module of any commercial FEA software like ABAQUS canbe helpful to generate the Points and subdomains in this proposed FPM. For example, Fig. 1(b)exhibits a typical point distribution and subdomain partition in FPM as converted form the ABAQUSmeshing.However, unlike in the FEM, the trial and test functions in the present FPM are point-based . Inthe primal FPM, in order to describe the high-order gradient-dependent behavior, the trial functionsfor the displacement and electric potential (( u h and φ h ) in each subdomain are written in the formof a third-order Taylor expansion at the corresponding internal point. For instance, in the FPMsubdomain E which contains an internal point P : u h ( x, y ) = u + ( x − x ) u , (cid:12)(cid:12)(cid:12) P + ( y − y ) u , (cid:12)(cid:12)(cid:12) P + 12 ( x − x ) u , (cid:12)(cid:12)(cid:12) P + ( x − x ) ( y − y ) u , (cid:12)(cid:12)(cid:12) P + 12 ( y − y ) u , (cid:12)(cid:12)(cid:12) P + 16 ( x − x ) u , (cid:12)(cid:12)(cid:12) P + 12 ( x − x ) ( y − y ) u , (cid:12)(cid:12)(cid:12) P + 12 ( x − x ) ( y − y ) u , (cid:12)(cid:12)(cid:12) P + 16 ( y − y ) u , (cid:12)(cid:12)(cid:12) P (26)13 a) (b)Figure 1: The domain Ω and its partitions. (a) Voronoi Diagram partition. (b) Quadrilateral and triangular partition(converted from FEA meshing). where [ x , y ] are the coordinates of point P . u is the value of u h at P . The first to third derivatives (cid:2) Du (cid:3) (cid:12)(cid:12)(cid:12) P = (cid:2) u T , , u T , , u T , , u T , , u T , , u T , , u T , , u T , , u T , (cid:3) T (cid:12)(cid:12)(cid:12) P are as yet unknown. Here weintroduce the local Differential Quadrature method to determine these high-order derivatives interms of the displacements at a few surrounding support Points. The conventional differential quadrature (DQ) method is a widely used numerical discretizationtechnique. However, usually based on 1D basis functions (e.g., polynomials), the conventional DQmethod can only be applied along a mesh-line and cannot be implemented directly with meshlessschemes. In recent years, several multi-dimensional DQ methods have been developed based onthe meshless Radial Basis Functions (RBFs). In the current work, we introduce and modify thelocal RBF-DQ method proposed by Shu et al. [52] in 2003 which approximates the 2D derivativesdirectly using a limited number of supporting points.In the current methodology, as shown in Fig. 1, for a given point P ∈ E (red), its supportis defined to involve all its nearest (blue) and second (green) neighboring points. Here the nearestneighboring points are the points in subdomains sharing boundaries with E , while the secondneighboring points are nearest neighbors of the former. For the points on or close to the domainboundary ∂ Ω , the third neighboring points are also taken into consideration to remedy the lack ofeffective supporting points. These supporting points are named as P , P , · · · , P m .14any RBFs can be used as basis function in the local RBF-DQ method. The multiquadric(MQ), inverse-MQ and Gaussians RBFs are all available to incorporate within the FPM approach,and achieve similar accuracy with appropriate parameters. Here we use the MQ-RBF as an example: ψ ( r ) = √ r + c , (27)where r is the radial length from the conference point. c = c d is a constant parameter, where d is the diameter of the minimal-diameter circle enclosing all the supporting Points.Thus, for supporting Points P i ( x i ) , the approximation of the displacement field u j ( x, y ) ( j = 1 , can be written as: u j ( x, y ) = m X i =0 λ i ψ ( k x − x i k ) + λ m +1 + λ m +2 x + λ m +3 y, (28)in which the coefficients λ i are determined by using the ( m + 1) collocating equations at P i and thefollowing conditions: m X i =0 λ i = 0 , m X i =0 λ i x i = 0 , m X i =0 λ i y i = 0 . (29)In the DQ method, any partial derivative of displacement at point P can be approximated by aweighted linear sum of the value of u h at all its supporting points: ∂ s + t ∂x s ∂y t u h (cid:12)(cid:12)(cid:12) P = m X i =0 W ( s,t ) i u h (cid:12)(cid:12) P i , (30)where W ( s,t ) i is the weighting coefficient corresponding to point P i . After rearrangement, the deriva-tives under study at point P can be approximated as: (cid:2) Du (cid:3) (cid:12)(cid:12)(cid:12) P = (cid:0) B ⊗ I × (cid:1) u E , (31)where u E = (cid:2) u , u , u , u , · · · , u m , u m (cid:3) T .u i = [ u i , u i ] T is the value of u h at P i . B is the matrix of weighting coefficients W ( s,t ) i , I × is a unitmatrix, and ⊗ denotes the Kronecker product. The weighting coefficient matrix B can be obtainedfrom: B = G − (cid:2) DG (cid:3) , (32)15here G = · · · x x · · · x m y y · · · y m g ( x , y ) g ( x , y ) · · · g ( x m , y m ) ... ... . . . ... g m ( x , y ) g m ( x , y ) · · · g m ( x m , y m ) , (cid:2) DG (cid:3) = g , ( x , y ) · · · g m, ( x , y )0 0 1 g , ( x , y ) · · · g m, ( x , y )0 0 0 g , ( x , y ) · · · g m, ( x , y )0 0 0 g , ( x , y ) · · · g m, ( x , y )0 0 0 g , ( x , y ) · · · g m, ( x , y )0 0 0 g , ( x , y ) · · · g m, ( x , y )0 0 0 g , ( x , y ) · · · g m, ( x , y )0 0 0 g , ( x , y ) · · · g m, ( x , y )0 0 0 g , ( x , y ) · · · g m, ( x , y ) T ,g i ( x, y ) = ψ i ( x, y ) − x y − x y + x i ( y − y ) + y i ( x − x )( x − x ) ( y − y ) − ( x − x ) ( y − y ) ψ ( x, y ) − x y − x y + x i ( y − y ) + y i ( x − x )( x − x ) ( y − y ) − ( x − x ) ( y − y ) ψ ( x, y ) − x y − x y + x i ( y − y ) + y i ( x − x )( x − x ) ( y − y ) − ( x − x ) ( y − y ) ψ ( x, y ) ,ψ i ( x, y ) = q ( x − x i ) + ( y − y i ) + c . Note that the accuracy of local RBF-DQ approximation is excellent for the first derivative;whereas, for higher-order approximations, the accuracy decreases. This implies that a mixed for-mulation may help to improve the accuracy of the current primal FPM.
Substituting the approximation of the derivatives into Eqn. 26, the trial function for displacement u h in the subdomain E is achieved (see Eqn. 33). Similarly, we can also generate the shape and16rial functions for the electric potential φ h : u h ( x, y ) = N u ( x, y ) u E , (33) φ h ( x, y ) = N φ ( x, y ) φ E , (34)where φ E = (cid:2) φ , φ , · · · , φ m (cid:3) T , N u ( x, y ) = N φ ( x, y ) ⊗ I × , N φ ( x, y ) = N ( x, y ) · B + [1 , , · · · , × ( m +1) , N ( x, y ) = h x − x , y − y ,
12 ( x − x ) , ( x − x ) ( y − y ) ,
12 ( y − y ) , · · ·
16 ( x − x ) ,
12 ( x − x ) ( y − y ) ,
12 ( x − x ) ( y − y ) ,
16 ( y − y ) i . The trial function and shape function are defined in each subdomain, following the same process.Since no continuity requirement is as yet enforced at the internal boundaries, the shape and trialfunctions can be discontinuous. Fig. 2(a) shows the graph of a shape function for the domain andpartition shown in Fig. 1(b) with 116 points. The trial function simulating an exponential function u a = e − [ ( x . +( y . ] / is presented in Fig. 2(b). As can be seen, the trial function is a cubicpolynomial in each subdomain. It is piecewise-continuous. However, as the shape function in eachsubdomain also depends on the values of the neighboring Points, the entire shape function showsa weakly continuous tendency. The test function for displacement v h and electrical potential τ h inthe Galerkin weak-form in the FPM are prescribed to possess the same shape as u h and φ h respec-tively. As a result of the discontinuous trial functions, the FPM also has a great natural potential inanalyzing systems involving cracks, ruptures and fragmentations.Unfortunately, if the above trial and test functions are applied in the traditional Galerkin weakform directly, the discontinuity will lead to inconsistent and inaccurate results. In order to resolvethat problem, we introduce the Numerical Flux Corrections to the present FPM.17 a) (b)Figure 2: The shape and trial function in the primal FPM. (a) The shape function. (b) The trial function for u a = e − [ ( x . +( y . ] / . We multiply the governing equations (Eqn. 3 and 4) by the test functions v hi and τ i respectivelyin each subdomain E , then apply the Gauss divergence theorem: Z E v i,j σ ij dΩ + Z E v i,jk µ ijk dΩ + Z E v i,j σ ESij dΩ= Z E v i b i dΩ + Z ∂E v i (cid:0) σ ij − µ ijk,k + σ ESij (cid:1) n j dΓ + Z ∂E v i,j µ ijk n k dΓ , (35) − Z E τ ,i D i dΩ − Z E τ ,ij Q ij dΩ (36) = Z E τ q dΩ − Z ∂E τ ( D i − Q ij,j ) n i dΓ − Z ∂E τ ,i Q ij n j dΓ . (37)where ∂E is the boundary of the subdomain E .Let Γ h denotes the set of all internal boundaries, and Γ = Γ h + ∂ Ω = Γ h + ∂ Ω u + ∂ Ω Q = Γ h + ∂ Ω φ + ∂ Ω ω = Γ h + ∂ Ω d + ∂ Ω R = Γ h + ∂ Ω P + ∂ Ω Z is the set of all internal and external boundaries.Each internal boundary e ∈ Γ h is shared by two subdomains E and E , i.e., e = ∂E ∩ ∂E (theorder does not affect the formulation of FPM). With n i being the unit vector normal to e and pointingoutward from E i , we define n e = n = − n for e ∈ Γ h . And for e ∈ ∂ Ω , n e = n , where n isthe outward unit normal vector of the only one neighboring subdomain E . By summing the above18quations over the entire domain, we rewrite them in the matrixâĂŹs forms: Z Ω ε T ( v ) σ ( u , φ ) dΩ + Z Ω κ T ( v ) µ ( u , φ ) dΩ + Z Ω ˆ ε T ( v ) σ ES ( u , φ ) dΩ − Z Γ h (cid:2)(cid:2) ˆ ε T ( v ) (cid:3)(cid:3) n e { µ ( u , φ ) } dΓ − Z Γ h (cid:8) ˆ ε T ( v ) (cid:9) n e [[ µ ( u , φ )]] dΓ − Z Γ h (cid:2)(cid:2) v T (cid:3)(cid:3) (cid:8) n e σ ( u , φ ) − n e µ , ( u , φ ) − n e µ , ( u , φ ) + n e σ ES ( u , φ ) (cid:9) dΓ − Z Γ h (cid:8) v T (cid:9) (cid:2)(cid:2) n e σ ( u , φ ) − n e µ , ( u , φ ) − n e µ , ( u , φ ) + n e σ ES ( u , φ ) (cid:3)(cid:3) dΓ= Z Ω v T b dΩ + Z ∂ Ω (cid:2)(cid:2) ˆ ε T ( v ) (cid:3)(cid:3) n e { µ ( u , φ ) } dΓ+ Z ∂ Ω (cid:2)(cid:2) v T (cid:3)(cid:3) (cid:8) n e σ ( u , φ ) − n e µ , ( u , φ ) − n e µ , ( u , φ ) + n e σ ES ( u , φ ) (cid:9) dΓ , (38) Z Ω E T ( τ ) D ( u , φ ) dΩ + Z Ω V T ( τ ) Q ( u , φ ) dΩ+ Z Γ h [[ τ ]] (cid:8) n e T D ( u , φ ) − n e Q , ( u , φ ) − n e Q , ( u , φ ) (cid:9) dΓ+ Z Γ h { τ } (cid:2)(cid:2) n e T D ( u , φ ) − n e Q , ( u , φ ) − n e Q , ( u , φ ) (cid:3)(cid:3) dΓ − Z Γ h (cid:2)(cid:2) E T ( τ ) (cid:3)(cid:3) n e { Q ( u , φ ) } dΓ − Z Γ h (cid:8) E T ( τ ) (cid:9) n e [[ Q ( u , φ )]] dΓ= Z Ω τ q dΩ − Z ∂ Ω [[ τ ]] (cid:8) n e T D ( u , φ ) − n e Q , ( u , φ ) − n e Q , ( u , φ ) (cid:9) dΓ+ Z ∂ Ω (cid:2)(cid:2) E T ( τ ) (cid:3)(cid:3) n e { Q ( u , φ ) } dΓ . (39)where ˆ ε T ( u ) = [ u , , u , , u , , u , ] T , (40)The jump operator [[ · ]] and average operator {·} are defined as (for ∀ w ): [[ w ]] = w (cid:12)(cid:12)(cid:12) E e − w (cid:12)(cid:12)(cid:12) E e e ∈ Γ h w (cid:12)(cid:12)(cid:12) e e ∈ ∂ Ω , { w } = (cid:18) w (cid:12)(cid:12)(cid:12) E e + w (cid:12)(cid:12)(cid:12) E e (cid:19) e ∈ Γ h w (cid:12)(cid:12)(cid:12) e e ∈ ∂ Ω . (41)The following boundary terms in Eqn. 38 and 39 can be divided into normal and tangent com-19onents and then integrated by parts: Z ∂ Ω (cid:2)(cid:2) ˆ ε T ( v ) (cid:3)(cid:3) n e { µ ( u , φ ) } dΓ = Z ∂ Ω d ∪ ∂ Ω R ˆ ε T ( v ) n e T4 n e n e µ ( u , φ ) dΓ − Z ∂ Ω u ∪ ∂ Ω Q v T c s e T4 s e n e µ , ( u , φ ) dΓ − Z ∂ Ω u ∪ ∂ Ω Q v T c s e T4 s e n e µ , ( u , φ ) dΓ+ X C h (cid:2)(cid:2) v T (cid:3)(cid:3) s e n e { µ ( u , φ ) } + X C h (cid:8) v T (cid:9) s e n e [[ µ ( u , φ )]] , (42) Z ∂ Ω (cid:2)(cid:2) E T ( τ ) (cid:3)(cid:3) n e { Q ( u , φ ) } dΓ = Z ∂ Ω P ∪ ∂ Ω Z E T ( τ ) n e n e T n e Q ( u , φ ) dΓ (43) + Z ∂ Ω φ ∪ ∂ Ω ω τ c s e s e T n e Q , ( u , φ ) dΓ + Z ∂ Ω φ ∪ ∂ Ω ω τ c s e s e T n e Q , ( u , φ ) dΓ (44) − X C h [[ τ ]] s e T n e { Q ( u , φ ) } − X C h { τ } s e T n e [[ Q ( u , φ )]] . (45)where C h denotes the set of all edges of the external boundaries, i.e., C h = S e ∈ ∂ Ω ∂e , s e is the unittangent vector pointing from E to E (the order is the same as the jump and average operatorsdefined in Eqn. 41), and s e · n e = 0 . The corresponding matrices in Eqn. 38 – 45 are shown inAppendix A.When ( u , φ ) are exact solutions, due to continuity conditions, we have: [[ µ ( u , φ )]] = , (cid:2)(cid:2) n e σ ( u , φ ) − n e µ , ( u , φ ) − n e µ , ( u , φ ) + n e σ ES ( u , φ ) (cid:3)(cid:3) = , (cid:2)(cid:2) n e T D ( u , φ ) − n e Q , ( u , φ ) − n e Q , ( u , φ ) (cid:3)(cid:3) = [[ Q ( u , φ )]] = , , for ∀ e ∈ Γ h . Hence, we can eliminate the corresponding terms in Eqn. 38 and 39. Similarly, [[ u ]] = , (cid:2)(cid:2) ˆ ε T ( u ) (cid:3)(cid:3) = , [[ φ ]] = 0 and (cid:2)(cid:2) E T ( φ ) (cid:3)(cid:3) = . We substitute Eqn. 15 into Eqn. 38. Then, by adding antithetic20erms and equivalent terms to both sides, the previous equations can be written in a symmetric form: Z Ω ε T ( v ) σ ( u , φ ) dΩ + Z Ω κ T ( v ) µ ( u , φ ) dΩ − Z Γ h ∪ ∂ Ω u (cid:2)(cid:2) v T (cid:3)(cid:3) n e σ ( u , φ ) − n e µ , ( u , φ ) − n e µ , ( u , φ )+ n e b D ( u , φ ) E ( φ ) + n e b Q ( u , φ ) V ( φ ) dΓ − Z Γ h ∪ ∂ Ω u (cid:8) σ T ( v , n e T1 − µ T , ( v , n e T21 − µ T , ( v , n e T22 (cid:9) [[ u ]] dΓ − Z Γ h ∪ ∂ Ω φ (cid:8) D T ( v , n e − Q T , ( v , n e T51 − Q T , ( v , n e T52 (cid:9) [[ φ ]] dΓ − Z Γ h (cid:2)(cid:2) ˆ ε T ( v ) (cid:3)(cid:3) n e { µ ( u , φ ) } dΓ − Z Γ h (cid:8) µ T ( v , (cid:9) n e T3 [[ ˆ ε ( u )]] dΓ+ Z Γ h (cid:8) Q T ( v , (cid:9) n e T6 [[ E ( φ )]] dΓ+ Z ∂ Ω u v T (cid:2) c s e T4 s e n e µ , ( u , φ ) + c s e T4 s e n e µ , ( u , φ ) (cid:3) dΓ+ Z ∂ Ω u (cid:2) µ T , ( v , n e T3 s e T4 s e c T1 + µ T , ( v , n e T3 s e T4 s e c T2 (cid:3) u dΓ+ Z ∂ Ω φ (cid:2) Q T , ( v , n e T6 s e s e T c T3 + Q T , ( v , n e T6 s e s e T c T4 (cid:3) φ dΓ − Z ∂ Ω d ˆ ε T ( v ) n e T4 n e n e µ ( u , φ ) dΓ − Z ∂ Ω d µ T ( v , n e T3 n e T4 n e ˆ ε ( u ) dΓ+ Z ∂ Ω P Q T ( v , n e T6 n e n e T E ( φ ) dΓ + X C h (cid:8) Q T ( v , (cid:9) n e T6 s e [[ φ ]] − X C h (cid:2)(cid:2) v T (cid:3)(cid:3) s e n e { µ ( u , φ ) } − X C h (cid:8) µ T ( v , (cid:9) n e T3 s e T4 [[ u ]]= Z Ω v T b dΩ − Z Ω ˆ ε T ( v ) h b D ( u , φ ) E ( φ ) + b Q ( u , φ ) V ( φ ) i dΩ+ Z ∂ Ω Q v T n e σ ( u , φ ) − (cid:0) n e + c s e T4 s e n e (cid:1) µ , ( u , φ ) − (cid:0) n e + c s e T4 s e n e (cid:1) µ , ( u , φ ) + n e σ ES ( u , φ ) dΓ+ Z ∂ Ω R ˆ ε T ( v ) n e T4 n e n e µ ( u , φ ) dΓ − Z ∂ Ω u σ T ( v , n e T1 − µ T , ( v , (cid:0) n e T21 + n e T3 s e T4 s e c T1 (cid:1) − µ T , ( v , (cid:0) n e T22 + n e T3 s e T4 s e c T2 (cid:1) u dΓ − Z ∂ Ω φ D T ( v , n e − Q T , ( v , (cid:0) n e T51 + n e T6 s e s e T c T3 (cid:1) − Q T , ( v , (cid:0) n e T52 + n e T6 s e s e T c T4 (cid:1) φ dΓ − Z ∂ Ω d µ T ( v , n e T3 n e T4 n e ˆ ε ( u ) dΓ + Z ∂ Ω P Q T ( v , n e T6 n e n e T E ( φ ) dΓ , (46)21 Z Ω E T ( τ ) D ( u , φ ) dΩ − Z Ω V T ( τ ) Q ( u , φ ) dΩ − Z Γ h ∪ ∂ Ω u σ T ( , τ ) n e T1 − µ T , ( , τ ) n e T21 − µ T , ( , τ ) n e T22 + E T ( τ ) b D T ( u , φ ) n e T4 + V T ( τ ) b Q T ( u , φ ) n e T4 [[ u ]] dΓ − Z Γ h ∪ ∂ Ω φ [[ τ ]] (cid:8) n e T D ( u , φ ) − n e Q , ( u , φ ) − n e Q , ( u , φ ) (cid:9) dΓ − Z Γ h ∪ ∂ Ω φ (cid:8) D T ( , τ ) n e − Q T , ( , τ ) n e T51 − Q T , ( , τ ) n e T52 (cid:9) [[ φ ]] dΓ+ Z Γ h (cid:2)(cid:2) E T ( τ ) (cid:3)(cid:3) n e { Q ( u , φ ) } dΓ + Z Γ h (cid:8) Q T ( , τ ) (cid:9) n e T6 [[ E ( φ )]] dΓ − Z Γ h (cid:8) µ T ( , τ ) (cid:9) n e T3 [[ ˆ ε ( u )]] dΓ+ Z ∂ Ω u (cid:2) µ T , ( , τ ) n e T3 s e T4 s e c T1 + µ T , ( , τ ) n e T3 s e T4 s e c T2 (cid:3) u dΓ+ Z ∂ Ω φ τ (cid:2) c s e s e T n e Q , ( u , φ ) + c s e s e T n e Q , ( u , φ ) (cid:3) dΓ+ Z ∂ Ω φ (cid:2) Q T , ( , τ ) n e T6 s e s e T c T3 + Q T , ( , τ ) n e T6 s e s e T c T4 (cid:3) φ dΓ+ Z ∂ Ω P E T ( τ ) n e n e T n e Q ( u , φ ) dΓ + Z ∂ Ω P Q T ( , τ ) n e T6 n e n e T E ( φ ) dΓ − Z ∂ Ω d µ T ( , τ ) n e T3 n e T4 n e ˆ ε ( u ) dΓ − X C h (cid:8) µ T ( , τ ) (cid:9) n e T3 s e T4 [[ u ]]+ X C h [[ τ ]] s e T n e { Q ( u , φ ) } + X C h (cid:8) Q T ( , τ ) (cid:9) n e T6 s e [[ φ ]]= − Z Ω τ q dΩ + Z ∂ Ω ω τ n e T D ( u , φ ) − (cid:0) n e + c s e s e T n e (cid:1) Q , ( u , φ ) − (cid:0) n e + c s e s e T n e (cid:1) Q , ( u , φ ) dΓ − Z ∂ Ω Z E T ( τ ) n e n e T n e Q ( u , φ ) dΓ − Z ∂ Ω u σ T ( , τ ) n e T1 − µ T , ( , τ ) (cid:0) n e T21 + n e T3 s e T4 s e c T1 (cid:1) − µ T , ( , τ ) (cid:0) n e T22 + n e T3 s e T4 s e c T2 (cid:1) + E T ( τ ) b D T ( u , φ ) n e T4 + V T ( τ ) b Q T ( u , φ ) n e T4 u dΓ − Z ∂ Ω φ D T ( , τ ) n e − Q T , ( , τ ) (cid:0) n e T51 + n e T6 s e s e T c T3 (cid:1) − Q T , ( , τ ) (cid:0) n e T52 + n e T6 s e s e T c T4 (cid:1) φ dΓ − Z ∂ Ω d µ T ( , τ ) n e T3 n e T4 n e ˆ ε ( u ) dΓ+ Z ∂ Ω P Q T ( , τ ) n e T6 n e n e T E ( φ ) dΓ . (47)22urthermore, in order to impose the essential boundary conditions ( ∂ Ω u , ∂ Ω d , ∂ Ω φ and ∂ Ω P )and the continuity condition across subdomains, we employ the Numerical Flux Correction, whichis widely used in Discontinuous Galerkin (DG) Methods [53, 54] to ensure the consistency andstability of the method. Here we use the Interior Penalty (IP) Numerical Flux Corrections. A set ofpenalty parameters are introduced, in which η = [ η , η , η , η ] and η = [ η , η , η , η ] arerelated to the essential boundary conditions and the interior continuity conditions respectively. TheFPM is only stable when the penalty parameters are large enough.Finally, after substituting the constitutive equations and the following boundary conditions (seeEqn. 6 – 13): u = e u , on ∂ Ω u , (48) φ = e φ, on ∂ Ω φ , (49) n e ˆ ε ( u ) = e d , on ∂ Ω d , (50) n e T E ( φ ) = − e P , on ∂ Ω P , (51) n e σ ( u , φ ) − (cid:0) n e + c s e T4 s e n e (cid:1) µ , ( u , φ ) − (cid:0) n e + c s e T4 s e n e (cid:1) µ , ( u , φ ) + n e σ ES ( u , φ ) = e Q , on ∂ Ω Q , (52) n e T D ( u , φ ) − (cid:0) n e + c s e s e T n e (cid:1) Q , ( u , φ ) − (cid:0) n e + c s e s e T n e (cid:1) Q , ( u , φ ) = − e ω, on ∂ Ω ω , (53) n e n e µ ( u , φ ) = e R , on ∂ Ω R , (54) n e T n e Q ( u , φ ) = e Z, on ∂ Ω Z , (55)23e can obtain the final symmetric formulation for the primal FPM for full flexoelectric theory: Z Ω ε T ( v ) C σε ε ( u ) dΩ + Z Ω κ T ( v ) C µκ κ ( u ) dΩ − Z Ω ε T ( v ) ( eE ( φ ) + bV ( φ )) dΩ − Z Ω κ T ( v ) aE ( φ ) dΩ − Z Γ h ∪ ∂ Ω u (cid:0)(cid:2)(cid:2) v T (cid:3)(cid:3) n e C σε { ε ( u ) } + (cid:8) ε T ( v ) (cid:9) C σε n e T1 [[ u ]] (cid:1) dΓ+ Z Γ h ∪ ∂ Ω u (cid:0)(cid:2)(cid:2) v T (cid:3)(cid:3) n e C µκ { κ , ( u ) } + (cid:8) κ T , ( v ) (cid:9) C µκ n e T21 [[ u ]] (cid:1) dΓ+ Z Γ h ∪ ∂ Ω u (cid:0)(cid:2)(cid:2) v T (cid:3)(cid:3) n e C µκ { κ , ( u ) } + (cid:8) κ T , ( v ) (cid:9) C µκ n e T22 [[ u ]] (cid:1) dΓ+ Z Γ h ∪ ∂ Ω u (cid:2)(cid:2) v T (cid:3)(cid:3) n(cid:16) n e b − n e ac − n e ac − n e b Q ( u , φ ) (cid:17) V ( φ ) o dΓ+ Z Γ h ∪ ∂ Ω u (cid:2)(cid:2) v T (cid:3)(cid:3) n(cid:16) n e e − n e b D ( u , φ ) (cid:17) E ( φ ) o dΓ − Z Γ h ∪ ∂ Ω φ (cid:8) ε T ( v ) (cid:9) en e [[ φ ]] dΓ − Z Γ h ∪ ∂ Ω φ (cid:8) κ T ( v ) (cid:9) (cid:0) an e − c T7 bn e T51 − c T8 bn e T52 (cid:1) [[ φ ]] dΓ − Z Γ h (cid:0)(cid:2)(cid:2) ˆ ε T ( v ) (cid:3)(cid:3) n e C µκ { κ ( u ) } + (cid:8) κ T ( v ) (cid:9) C µκ n e T3 [[ ˆ ε ( u )]] (cid:1) dΓ+ Z Γ h (cid:2)(cid:2) ˆ ε T ( v ) (cid:3)(cid:3) n e a { E ( φ ) } dΓ + Z Γ h (cid:8) ε T ( v ) (cid:9) bn e T6 [[ E ( φ )]] dΓ+ Z ∂ Ω u (cid:0) v T c s e T4 s e n e C µκ κ , ( u ) + κ T , ( v ) C µκ n e T3 s e T4 s e c T1 u (cid:1) dΓ+ Z ∂ Ω u (cid:0) v T c s e T4 s e n e C µκ κ , ( u ) + κ T , ( v ) C µκ n e T3 s e T4 s e c T2 u (cid:1) dΓ − Z ∂ Ω u v T (cid:0) c s e T4 s e n e ac + c s e T4 s e n e ac (cid:1) V ( φ ) dΓ+ Z ∂ Ω φ κ T ( v ) (cid:0) c T7 bn e T6 s e s e T c T3 + c T8 bn e T6 s e s e T c T4 (cid:1) φ dΓ − Z ∂ Ω d (cid:0) ˆ ε T ( v ) n e T4 n e n e C µκ κ ( u ) + κ T ( v ) C µκ n e T3 n e T4 n e ˆ ε ( u ) (cid:1) dΓ+ Z ∂ Ω d ˆ ε T ( v ) n e T4 n e n e aE ( φ ) dΓ + Z ∂ Ω P ε T ( v ) bn e T6 n e n e T E ( φ ) dΓ+ Z ∂ Ω u η h e v T u dΓ + Z Γ h η h e (cid:2)(cid:2) v T (cid:3)(cid:3) [[ u ]] dΓ+ Z ∂ Ω d η h e ˆ ε T ( v ) n e T4 n e ˆ ε ( u ) dΓ + Z Γ h η h e (cid:2)(cid:2) ˆ ε T ( v ) (cid:3)(cid:3) n e T4 n e [[ ˆ ε ( u )]] dΓ (56)24 X C h (cid:2)(cid:2) v T (cid:3)(cid:3) s e n e C µκ { κ ( u ) } − X C h (cid:8) κ T ( v ) (cid:9) C µκ n e T3 s e T4 [[ u ]]+ X C h (cid:8) ε T ( v ) (cid:9) bn e T6 s e [[ φ ]] + X C h (cid:2)(cid:2) v T (cid:3)(cid:3) s e n e a { E ( φ ) } = Z Ω v T b dΩ − Z Ω ˆ ε T ( v ) (cid:16) b D ( u , φ ) E ( φ ) + b Q ( u , φ ) V ( φ ) (cid:17) dΩ+ Z ∂ Ω Q v T e Q dΓ + Z ∂ Ω R ˆ ε T ( v ) n e T4 e R dΓ − Z ∂ Ω u ε T ( v ) C σε n e T1 − κ T , ( v ) C µκ (cid:0) n e T21 + n e T3 s e T4 s e c T1 (cid:1) − κ T , ( v ) C µκ (cid:0) n e T22 + n e T3 s e T4 s e c T2 (cid:1) e u dΓ − Z ∂ Ω φ ε T ( v ) en e + κ T ( v ) (cid:2) an e − c T7 b (cid:0) n e T51 + n e T6 s e s e T c T3 (cid:1) − c T8 b (cid:0) n e T52 + n e T6 s e s e T c T4 (cid:1)(cid:3) e φ dΓ − Z ∂ Ω d κ T ( v ) C µκ n e T3 n e T4 e d dΓ − Z ∂ Ω P ε T ( v ) bn e T6 n e e P dΓ+ Z ∂ Ω u η h e v T e u dΓ + Z ∂ Ω d η h e ˆ ε T ( v ) n e T4 e d dΓ , Z Ω E T ( τ ) ΛE ( φ ) dΩ − Z Ω V T ( τ ) ΦV ( φ ) dΩ − Z Ω (cid:0) E T ( τ ) e T + V T ( τ ) b T (cid:1) ε ( u ) dΩ − Z Ω E T ( τ ) a T κ ( u ) dΩ+ Z Γ h ∪ ∂ Ω u n E T ( τ ) (cid:16) e T n e T1 − b D T ( u , φ ) n e T4 (cid:17)o [[ u ]] dΓ+ Z Γ h ∪ ∂ Ω u n V T ( τ ) (cid:16) b T n e T1 − c T5 a T n e T21 − c T6 a T n e T22 − b Q T ( u , φ ) n e T4 (cid:17)o [[ u ]] dΓ − Z Γ h ∪ ∂ Ω φ (cid:0) [[ τ ]] n e T Λ { E ( φ ) } + (cid:8) E T ( τ ) (cid:9) Λn e [[ φ ]] (cid:1) dΓ+ Z Γ h ∪ ∂ Ω φ (cid:0) [[ τ ]] n e Φ { V , ( φ ) } + (cid:8) V T , ( τ ) (cid:9) Φn e T51 [[ φ ]] (cid:1) dΓ+ Z Γ h ∪ ∂ Ω φ (cid:0) [[ τ ]] n e Φ { V , ( φ ) } + (cid:8) V T , ( τ ) (cid:9) Φn e T52 [[ φ ]] (cid:1) dΓ − Z Γ h ∪ ∂ Ω φ [[ τ ]] n e T e T { ε ( u ) } dΓ − Z Γ h ∪ ∂ Ω φ [[ τ ]] (cid:0) n e T a T − n e b T c − n e b T c (cid:1) { κ ( u ) } dΓ+ Z Γ h (cid:0)(cid:2)(cid:2) E T ( τ ) (cid:3)(cid:3) n e Φ { V ( φ ) } + (cid:8) V T ( τ ) (cid:9) Φn e T6 [[ E ( φ )]] (cid:1) dΓ+ Z Γ h (cid:2)(cid:2) E T ( τ ) (cid:3)(cid:3) n e b T { ε ( u ) } dΓ + Z Γ h (cid:8) E T ( τ ) (cid:9) a T n e T3 [[ ˆ ε ( u )]] dΓ − Z ∂ Ω u V T ( τ ) (cid:0) c T5 a T n e T3 s e T4 s e c T1 + c T6 a T n e T3 s e T4 s e c T2 (cid:1) u dΓ+ Z ∂ Ω φ (cid:0) τ c s e s e T n e ΦV , ( φ ) + V T , ( τ ) Φn e T6 s e s e T c T3 φ (cid:1) dΓ+ Z ∂ Ω φ (cid:0) τ c s e s e T n e ΦV , ( φ ) + V T , ( τ ) Φn e T6 s e s e T c T4 φ (cid:1) dΓ+ Z ∂ Ω φ τ (cid:0) c s e s e T n e b T c + c s e s e T n e b T c (cid:1) κ ( u ) dΓ+ Z ∂ Ω d E T ( τ ) a T n e T3 n e T4 n e ˆ ε ( u ) dΓ + Z ∂ Ω P E T ( τ ) n e n e T n e b T ε ( u ) dΓ+ Z ∂ Ω P (cid:0) E T ( τ ) n e n e T n e ΦV ( φ ) + V T ( τ ) Φn e T6 n e n e T E ( φ ) (cid:1) dΓ+ Z ∂ Ω φ η h e τ φ dΓ + Z Γ h η h e [[ τ ]] [[ φ ]] dΓ+ Z ∂ Ω P η h e E T ( τ ) n e n e T E ( φ ) dΓ + Z Γ h η h e (cid:2)(cid:2) E T ( τ ) n e (cid:3)(cid:3) n e T4 n e (cid:2)(cid:2) n e T E ( φ ) (cid:3)(cid:3) dΓ (57)26 X C h [[ τ ]] s e T n e Φ { V ( φ ) } + X C h (cid:8) V T ( τ ) (cid:9) Φn e T6 s e [[ φ ]]+ X C h [[ τ ]] s e T n e b T { ε ( u ) } + X C h (cid:8) E T ( τ ) (cid:9) a T n e T3 s e T4 [[ u ]]= − Z Ω τ q dΩ − Z ∂ Ω ω τ e ω dΓ − Z ∂ Ω Z E T ( τ ) n e e Z dΓ+ Z ∂ Ω u E T ( τ ) (cid:16) e T n e T1 − b D T ( u , φ ) n e T4 (cid:17) e u dΓ+ Z ∂ Ω u V T ( τ ) b T n e T1 − c T5 a T (cid:0) n e T21 + n e T3 s e T4 s e c T1 (cid:1) − c T6 a T (cid:0) n e T22 + n e T3 s e T4 s e c T2 (cid:1) − b Q T ( u , φ ) n e T4 e u dΓ − Z ∂ Ω φ E T ( τ ) Λn e − V T , ( τ ) Φ (cid:0) n e T51 + n e T6 s e s e T c T3 (cid:1) − V T , ( τ ) Φ (cid:0) n e T52 + n e T6 s e s e T c T4 (cid:1) e φ dΓ+ Z ∂ Ω d E T ( τ ) a T n e T3 n e T4 e d dΓ − Z ∂ Ω P V T ( τ ) Φn e T6 n e e P dΓ+ Z ∂ Ω φ η h e τ e φ dΓ − Z ∂ Ω P η h e E T ( τ ) n e e P dΓ . where the normal and constant matrices n i and c i can be seen in Appendix A. h e is a boundary-dependent parameter with the unit of length. Here it is defined as the distance between the inter-nal points in subdomains sharing the boundary for e ∈ Γ h , and the smallest distance between thecentroid of the subdomain and the external boundary for e ∈ ∂ Ω . The penalty parameters are in-dependent of the boundary size, in which ( η , η , η , η ) have the same unit as of the Young’smodulus, and ( η , η , η , η ) have the same unit as 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 element shape functions are completely independent,a small penalty parameter η is enough to stabilize the present FPM. On the other hand, the accuracymay decrease if η is too large. η on the other hand, should be larger, which corresponds to theessential boundary conditions.It should be noted that the essential boundary conditions can also be imposed in an alternativeway in the FPM. When boundary points are employed, i.e., for ∂E i ∩ ∂ Ω = ∅ , P i ∈ ( ∂E i ∩ ∂ Ω) ,we can enforce u = e u and φ = e φ strongly at the boundary points and thus the Internal Penalty(IP) terms relating to η and η in Eqn. 56 and 57 vanish. However, when high-order essentialboundary conditions are under study ( ∂ Ω d ∪ ∂ Ω P = ∅ ), the IP terms are still required. Note that27he IP Numerical Flux Corrections can enforce the essential boundary conditions with any pointdistributions, regardless of whether there are points distributed on the boundary.At last, the formula of the FPM can be converted into the following form: K Puu (cid:0) u , φ (cid:1) K Puφ (cid:0) u , φ (cid:1)(cid:0) K Puφ (cid:0) u , φ (cid:1)(cid:1) T K Pφφ (cid:0) u , φ (cid:1) u φ = f Pu (cid:0) u , φ (cid:1) f Pφ (cid:0) u , φ (cid:1) , or K P ( x ) x = f P ( x ) , (58)where (cid:0) u , φ (cid:1) are the global nodal displacement and electric potential vectors. K P is the global stiff-ness matrix for the primal FPM, assembled by a series of point stiffness matrices ( K PE ) and boundarystiffness matrices ( K Ph , K Pu , K PQ , K Pd , K PR , K Pφ , K Pω , K PP and K PZ corresponding to Γ h , ∂ Ω u , ∂ Ω Q , ∂ Ω d , ∂ Ω R , ∂ Ω φ , ∂ Ω ω , ∂ Ω P and ∂ Ω Z respectively). f P is the global load vector. When evaluat-ing these stiffness matrices and load vectors, we use simple Gaussian quadrature rules. Numericalexamples have shown that only one integration point in each subdomain and only one integrationpoint on each internal and external boundary can result in good solutions with acceptable accura-cies. However, × or more integration points may help to further increase the accuracy, but onlyone point integration schemes are employed in the present paper. Further numerical implementationin details are omitted here.It should be pointed out that Eqn. 58 based on the full flexoelectric theory is nonlinear, i.e.,the global stiffness matrix K P and load vector f P are functions of the unknown displacement andelectric potential fields. As can be seen from Eqn. 56 and 57, the Point and boundary stiffnessmatrices and load vectors are dependent on e D ( u , φ ) and e Q ( u , φ ) in each subdomain. As a result,the global stiffness matrix and load vector has a linear relationship with x . The Newton-Raphsonmethod is then applied to solve the nonlinear equation. Notice that when the initial guess x (0) = , the first approximation x (1) = (cid:2) K P ( ) (cid:3) − f P ( ) is equivalent to the linear solution when theelectrostatic stress σ ES is neglected. In the reduced formulation (Eqn. 16 – 23), both the generalized electrostatic stress σ ES andelectric field gradient V are neglected. Thus the trial and test functions of the electric potential canbe written in a quadratic form: φ h ( x, y ) = N φ ( x, y ) φ E , N φ ( x, y ) = N R ( x, y ) · B + [1 , , · · · , × ( m +1) , N R ( x, y ) = (cid:20) x − x , y − y ,
12 ( x − x ) , ( x − x ) ( y − y ) ,
12 ( y − y ) , , , , (cid:21) . (59)28ollowing the same steps in section 3.2.1, we can obtain the linear primal FPM formula for thereduced flexoelectric theory: Z Ω ε T ( v ) C σε ε ( u ) dΩ + Z Ω κ T ( v ) C µκ κ ( u ) dΩ − Z Ω (cid:0) ε T ( v ) e + κ T ( v ) a (cid:1) E ( φ ) dΩ − Z Γ h ∪ ∂ Ω u (cid:0)(cid:2)(cid:2) v T (cid:3)(cid:3) n e C σε { ε ( u ) } + (cid:8) ε T ( v ) (cid:9) C σε n e T1 [[ u ]] (cid:1) dΓ+ Z Γ h ∪ ∂ Ω u (cid:0)(cid:2)(cid:2) v T (cid:3)(cid:3) n e C µκ { κ , ( u ) } + (cid:8) κ T , ( v ) (cid:9) C µκ n e T21 [[ u ]] (cid:1) dΓ+ Z Γ h ∪ ∂ Ω u (cid:0)(cid:2)(cid:2) v T (cid:3)(cid:3) n e C µκ { κ , ( u ) } + (cid:8) κ T , ( v ) (cid:9) C µκ n e T22 [[ u ]] (cid:1) dΓ − Z Γ h ∪ ∂ Ω u (cid:2)(cid:2) v T (cid:3)(cid:3) { n e aE , ( φ ) + n e aE , ( φ ) } dΓ+ Z Γ h ∪ ∂ Ω u (cid:2)(cid:2) v T (cid:3)(cid:3) n e e { E ( φ ) } dΓ + Z Γ h (cid:2)(cid:2) ˆ ε T ( v ) (cid:3)(cid:3) n e a { E ( φ ) } dΓ − Z Γ h (cid:0)(cid:2)(cid:2) ˆ ε T ( v ) (cid:3)(cid:3) n e C µκ { κ ( u ) } + (cid:8) κ T ( v ) (cid:9) C µκ n e T3 [[ ˆ ε ( u )]] (cid:1) dΓ − Z Γ h ∪ ∂ Ω φ (cid:8) ε T ( v ) (cid:9) en e [[ φ ]] dΓ − Z Γ h ∪ ∂ Ω φ (cid:8) κ T ( v ) (cid:9) an e [[ φ ]] dΓ+ Z ∂ Ω u (cid:0) v T c s e T4 s e n e C µκ κ , ( u ) + κ T , ( v ) C µκ n e T3 s e T4 s e c T1 u (cid:1) dΓ+ Z ∂ Ω u (cid:0) v T c s e T4 s e n e C µκ κ , ( u ) + κ T , ( v ) C µκ n e T3 s e T4 s e c T2 u (cid:1) dΓ − Z ∂ Ω u v T (cid:0) c s e T4 s e n e aE , ( φ ) + c s e T4 s e n e aE , ( φ ) (cid:1) dΓ − Z ∂ Ω d (cid:0) ˆ ε T ( v ) n e T4 n e n e C µκ κ ( u ) + κ T ( v ) C µκ n e T3 n e T4 n e ˆ ε ( u ) (cid:1) dΓ+ Z ∂ Ω d ˆ ε T ( v ) n e T4 n e n e aE ( φ ) dΓ + Z ∂ Ω u η h e v T u dΓ + Z Γ h η h e (cid:2)(cid:2) v T (cid:3)(cid:3) [[ u ]] dΓ+ Z ∂ Ω d η h e ˆ ε T ( v ) n e T4 n e ˆ ε ( u ) dΓ + Z Γ h η h e (cid:2)(cid:2) ˆ ε T ( v ) (cid:3)(cid:3) n e T4 n e [[ ˆ ε ( u )]] dΓ − X C h (cid:2)(cid:2) v T (cid:3)(cid:3) s e n e C µκ { κ ( u ) } − X C h (cid:8) κ T ( v ) (cid:9) C µκ n e T3 s e T4 [[ u ]]+ X C h (cid:2)(cid:2) v T (cid:3)(cid:3) s e n e a { E ( φ ) } (60)29 Z Ω v T b dΩ + Z ∂ Ω Q v T e Q dΓ + Z ∂ Ω R ˆ ε T ( v ) n e T4 e R dΓ − Z ∂ Ω u ε T ( v ) C σε n e T1 − κ T , ( v ) C µκ (cid:0) n e T21 + n e T3 s e T4 s e c T1 (cid:1) − κ T , ( v ) C µκ (cid:0) n e T22 + n e T3 s e T4 s e c T2 (cid:1) e u dΓ − Z ∂ Ω φ (cid:0) ε T ( v ) en e + κ T ( v ) an e (cid:1) e φ dΓ − Z ∂ Ω d κ T ( v ) C µκ n e T3 n e T4 e d dΓ+ Z ∂ Ω u η h e v T e u dΓ + Z ∂ Ω d η h e ˆ ε T ( v ) n e T4 e d dΓ , − Z Ω E T ( τ ) ΛE ( φ ) dΩ − Z Ω E T ( τ ) (cid:0) e T ε ( u ) + a T κ ( u ) (cid:1) dΩ − Z Γ h ∪ ∂ Ω u (cid:8) E T , ( τ ) a T n e T21 + E T , ( τ ) a T n e T22 (cid:9) [[ u ]] dΓ+ Z Γ h ∪ ∂ Ω u (cid:8) E T ( τ ) (cid:9) e T n e T1 [[ u ]] dΓ + Z Γ h (cid:8) E T ( τ ) (cid:9) a T n e T3 [[ ˆ ε ( u )]] dΓ − Z Γ h ∪ ∂ Ω φ [[ τ ]] n e T e T { ε ( u ) } dΓ − Z Γ h ∪ ∂ Ω φ [[ τ ]] n e T a T { κ ( u ) } dΓ − Z Γ h ∪ ∂ Ω φ (cid:0) [[ τ ]] n e T Λ { E ( φ ) } + (cid:8) E T ( τ ) (cid:9) Λn e [[ φ ]] (cid:1) dΓ − Z ∂ Ω u (cid:0) E T , ( τ ) a T n e T3 s e T4 s e c T1 + E T , ( τ ) a T n e T3 s e T4 s e c T2 (cid:1) u dΓ+ Z ∂ Ω d E T ( τ ) a T n e T3 n e T4 n e ˆ ε ( u ) dΓ + Z ∂ Ω φ η h e τ φ dΓ + Z Γ h η h e [[ τ ]] [[ φ ]] dΓ+ X C h (cid:8) E T ( τ ) (cid:9) a T n e T3 s e T4 [[ u ]]= − Z Ω τ q dΩ − Z ∂ Ω ω τ e ω dΓ + Z ∂ Ω u E T ( τ ) e T n e T1 e u dΓ − Z ∂ Ω u (cid:2) E T , ( τ ) a T (cid:0) n e T21 + n e T3 s e T4 s e c T1 (cid:1) + E T , ( τ ) a T (cid:0) n e T22 + n e T3 s e T4 s e c T2 (cid:1)(cid:3) e u dΓ − Z ∂ Ω φ E T ( τ ) Λn e e φ dΓ + Z ∂ Ω d E T ( τ ) a T n e T3 n e T4 e d dΓ + Z ∂ Ω φ η h e τ e φ dΓ . (61)In a matrix form: K Puu K Puφ (cid:0) K Puφ (cid:1) T K Pφφ u φ = f Pu f Pφ , or K P x = f P . (62)In the reduced theory, the global stiffness matrix and load vector K P and f P are independent of thedisplacement and electric potential field, i.e., Eqn. 62 is linear. Besides, the final stiffness matrix is30parse and symmetric, which is beneficial for modeling complex problems with a large number ofDoFs.Whereas the trial and test functions for displacement given in this section are cubic polynomial,numerical study shows that simplified quadratic trial and test functions can also by employed andcan achieve good accuracy (with relative error < κ ,i ( i = 1 , on the internal boundaries in Eqn. 56 – 61are neglected. Similarly, the trial and test functions for the electric potential can also be simplifiedto quadratic (in the full theory) or linear (in the reduced theory). And the corresponding termscontaining V ,i ( i = 1 , or E ,i ( i = 1 , in the above formulations are then absent. The cubic trialand test functions can improve the accuracy of the primal FPM with more terms in the TaylorâĂŹsexpansion.
4. Mixed Fragile Points Method (mixed FPM)
In this section, we develop a mixed FPM in which the displacement u , electric potential φ ,high-order stress µ and high-order electric displacement Q are all interpolated independently.Here we introduce the following independent variables and their trial functions in subdomain E : u h ( x, y ) = u + ( x − x ) u , (cid:12)(cid:12) P + ( y − y ) u , (cid:12)(cid:12) P = N Mu ( x, y ) u E , (63) φ h ( x, y ) = φ + ( x − x ) φ , (cid:12)(cid:12) P + ( y − y ) φ , (cid:12)(cid:12) P = N Mφ ( x, y ) φ E , (64) µ h ( x, y ) = µ + ( x − x ) µ , (cid:12)(cid:12) P + ( y − y ) µ , (cid:12)(cid:12) P = N Mµ ( x, y ) µ E , (65) Q h ( x, y ) = Q + ( x − x ) Q , (cid:12)(cid:12) P + ( y − y ) Q , (cid:12)(cid:12) P = N MQ ( x, y ) Q E . (66)where µ E = (cid:2) µ , µ , · · · , µ m T (cid:3) T , Q E = (cid:2) Q , Q , · · · , Q m T (cid:3) T , N Mu ( x, y ) = N Mφ ( x, y ) ⊗ I × , N Mµ ( x, y ) = N Mφ ( x, y ) ⊗ I × , N MQ ( x, y ) = N Mφ ( x, y ) ⊗ I × , N Mφ ( x, y ) = N M ( x, y ) · B M + [1 , , · · · , × ( m +1) , N M ( x, y ) = [ x − x , y − y ] , B M = B · [ I × , × ] T . (67)31imilar to the primal FPM, the trial functions in the mixed FPM are also piecewise-continuous.Yet for simplicity, these trial functions are just written as a linear Taylor expansion at the corre-sponding Fragile Point and all the high-order terms are omitted. The graph of the linear piecewise-continuous shape function is shown in Fig. 3(a). And the trial function simulating an exponentialfunction u a = e − √ ( x . +( y . is presented in Fig. 3(b). There are various approaches to ap-proximate the first derivatives based on the values of the supporting points. In addition to the localRBF-DQ method shown in section 3.1.2, Generalized Finite Difference (GFD) method can also beapplied to generate B M in Eqn. 67. The corresponding matrices can be seen in [40] and are omit-ted here. The test functions for displacement v h , electric potential τ h , high-order stress w h andhigh-order electric displacement ξ h have the same shape as their corresponding trial functions inthe mixed FPM. (a) (b)Figure 3: The shape and trial function in the mixed FPM. (a) The shape function. (b) The trial function for u a = e − √ ( x . +( y . . .2. Weak-form formulation4.2.1. Full theory First, we multiply the governing equations (Eqn. 3 and 4) by the test functions v hi and τ i respec-tively and integrate in subdomain E by parts: Z E ε T ( v ) σ ( u , φ, Q ) dΩ + Z E ˆ ε T ( v ) σ ES ( u , φ, Q ) dΩ − Z E ˆ ε T ( v ) c T9 µ , dΩ − Z E ˆ ε T ( v ) c T10 µ , dΩ= Z E v T b dΩ + Z ∂E v T (cid:0) n e σ ( u , φ, Q ) − n e µ , − n e µ , + n e σ ES ( u , φ, Q ) (cid:1) dΓ , (68) Z E E T ( τ ) D ( u , φ, µ ) dΩ − Z E E T ( τ ) c Q , dΩ − Z E E T ( τ ) c Q , dΩ= Z E τ q dΩ − Z ∂E τ (cid:0) n e T D ( u , φ, µ ) − n e Q , − n e Q , (cid:1) dΓ . (69)Similarly, the constitutive relations for µ and Q (see Eqn. 14) are also multiplied by their corre-sponding test functions w and ξ and then integrated by parts in E : Z E w T C − µκ µ dΩ + Z E w T , c ˆ ε ( u ) dΩ + Z E w T , c ˆ ε ( u ) dΩ − Z E w T C − µκ aE ( φ ) dΩ = Z ∂E w T n e T3 ˆ ε ( u ) dΓ , (70) Z E ξ T Φ − Q dΩ + Z E ξ T , c T5 E ( φ ) dΩ + Z E ξ T , c T6 E ( φ ) dΩ − Z E ξ T Φ − b T ε ( u ) dΩ = Z ∂E ξ T n e T6 E ( φ ) dΓ . (71)The high-order boundary conditions for µ and Q on ∂ Ω R and ∂ Ω Z are enforced by Lagrangemultipliers λ R = [ λ R , λ R ] and λ Z respectively. These Lagrange multipliers are defined on ∂ Ω R and ∂ Ω Z with constant shape function in each subdomain boundary: λ hR ( x, y ) = λ iR , λ hZ ( x, y ) = λ iZ , for ( x, y ) ∈ e i . (72)The nodal Lagrange multiplier vectors can be written as λ R = [ λ R , λ R , λ R , λ R , · · · , λ n R R , λ n R R ] T and λ Z = [ λ Z , λ Z , · · · , λ n Z Z ] T , where n R and n Z are the number of subdomain boundaries on ∂ Ω R and ∂ Ω Z respectively. The test functions of the Lagrange multipliers ˆ λ R and ˆ λ Z take the same shapeas λ R and λ Z . 33ollowing the same steps in the primal FPM, we can obtain the basic symmetric formula for themixed FPM: Z Ω ε T ( v ) (cid:0) C σε + bΦ − b T (cid:1) ε ( u ) dΩ − Z Ω ε T ( v ) eE ( φ ) dΩ − Z Ω ε T ( v ) bΦ − Q dΩ − Z Ω ˆ ε T ( v ) c T9 µ , dΩ − Z Ω ˆ ε T ( v ) c T10 µ , dΩ − Z Γ h ∪ ∂ Ω u (cid:2)(cid:2) v T (cid:3)(cid:3) n(cid:16) n e (cid:0) C σε + bΦ − b T (cid:1) − n e b QΦ − b T (cid:17) ε ( u ) o dΓ − Z Γ h ∪ ∂ Ω u n ε T ( v ) (cid:16)(cid:0) C σε + bΦ − b T (cid:1) n e T1 − bΦ − b Q T n e T4 (cid:17)o [[ u ]] dΓ+ Z Γ h ∪ ∂ Ω u (cid:2)(cid:2) v T (cid:3)(cid:3) n(cid:16) n e e − n e b D ( u , φ, µ ) (cid:17) E ( φ ) o dΓ+ Z Γ h ∪ ∂ Ω u (cid:2)(cid:2) v T (cid:3)(cid:3) { n e µ , + n e µ , } dΓ + Z Γ h (cid:8) ˆ ε T ( v ) (cid:9) n e [[ µ ]] dΓ+ Z Γ h ∪ ∂ Ω u (cid:2)(cid:2) v T (cid:3)(cid:3) n(cid:16) n e b − n e b Q (cid:17) Φ − Q o dΓ − Z Γ h ∪ ∂ Ω φ (cid:8) ε T ( v ) (cid:9) en e [[ φ ]] dΓ − Z ∂ Ω Q v T (cid:0) c s e T4 s e µ , + c s e T4 s e µ , (cid:1) dΓ + Z ∂ Ω R ˆ ε T ( v ) n e T4 n e n e µ dΓ+ Z ∂ Ω u η h e v T u dΓ + Z Γ h η h e (cid:2)(cid:2) v T (cid:3)(cid:3) [[ u ]] dΓ + X C h (cid:8) v T (cid:9) s e n e [[ µ ]] (73) = Z Ω v T b dΩ + Z ∂ Ω Q v T e Q dΓ − Z Ω ˆ ε T ( v ) b D ( u , φ, µ ) E ( φ ) dΩ − Z Ω ˆ ε T ( v ) b QΦ − Q dΩ+ Z Ω ˆ ε T ( v ) b QΦ − b T ε ( u ) dΩ − Z ∂ Ω φ ε T ( v ) en e e φ dΓ + Z ∂ Ω R ˆ ε T ( v ) n e T4 e R dΓ − Z ∂ Ω u ε T ( v ) (cid:16)(cid:0) C σε + bΦ − b T (cid:1) n e T1 − bΦ − b Q T n e T4 (cid:17) e u dΓ + Z ∂ Ω u η h e v T e u dΓ , Z Ω E T ( τ ) (cid:0) Λ + a T C − µκ a (cid:1) E ( φ ) dΩ − Z Ω E T ( τ ) e T ε ( u ) dΩ − Z Ω E T ( τ ) a T C − µκ µ dΩ + Z Ω E T ( τ ) n Q , dΩ + Z Ω E T ( τ ) n Q , dΩ − Z Γ h ∪ ∂ Ω φ [[ τ ]] n e T (cid:0) Λ + a T C − µκ a (cid:1) { E ( φ ) } dΓ − Z Γ h ∪ ∂ Ω φ [[ τ ]] n e T e T { ε ( u ) } dΓ − Z Γ h ∪ ∂ Ω φ (cid:8) E T ( τ ) (cid:9) (cid:0) Λ + a T C − µκ a (cid:1) n e [[ φ ]] dΓ − Z Γ h ∪ ∂ Ω φ [[ τ ]] n e T a T C − µκ { µ } dΓ+ Z Γ h ∪ ∂ Ω φ [[ τ ]] n e { Q , } dΓ + Z Γ h ∪ ∂ Ω φ [[ τ ]] n e { Q , } dΓ+ Z Γ h ∪ ∂ Ω u n E T ( τ ) (cid:16) e T n e T1 − b D T ( u , φ, µ ) n e T4 (cid:17)o [[ u ]] dΓ − Z Γ h (cid:8) E T ( τ ) (cid:9) n e [[ Q ]] dΓ + Z Γ h η h e [[ τ ]] [[ φ ]] dΓ − Z ∂ Ω Z E T ( τ ) n e n e T n e Q dΓ − Z ∂ Ω ω τ (cid:0) c s e s e T n e Q , + c s e s e T n e Q , (cid:1) dΓ + Z ∂ Ω φ η h e τ φ dΓ + X C h { τ } s e T n e [[ Q ]] (74) = − Z Ω τ q dΩ − Z ∂ Ω ω τ e ω dΓ + Z ∂ Ω u E T ( τ ) (cid:16) e T n e T1 − b D T ( u , φ, µ ) n e T4 (cid:17) e u dΓ − Z ∂ Ω Z E T ( τ ) n e e Z dΓ − Z ∂ Ω φ E T ( τ ) (cid:0) Λ + a T C − µκ a (cid:1) n e e φ dΓ + Z ∂ Ω φ η h e τ e φ dΓ , − Z Ω w T C − µκ µ dΩ − Z Ω w T C − µκ aE ( φ ) dΩ − Z Ω w T , c ˆ ε ( u ) dΩ − Z Ω w T , c ˆ ε ( u ) dΩ+ Z Γ h ∪ ∂ Ω u (cid:8) w T , n e T21 + w T , n e T22 (cid:9) [[ u ]] dΓ − Z Γ h ∪ ∂ Ω φ (cid:8) w T (cid:9) C − µκ an e [[ φ ]] dΓ+ Z Γ h (cid:2)(cid:2) w T (cid:3)(cid:3) n e T3 { ˆ ε ( u ) } dΓ + Z ∂ Ω R w T n e T3 n e T4 n e ˆ ε ( u ) dΓ + Z ∂ Ω R w T n e T3 n e T4 λ R dΓ − Z ∂ Ω Q w T , n e T3 s e T4 s e c T1 u dΓ − Z ∂ Ω Q w T , n e T3 s e T4 s e c T2 u dΓ + X C h (cid:2)(cid:2) w T (cid:3)(cid:3) n e T3 s e T4 { u } (75) = − Z ∂ Ω d w T n e T3 n e T4 e d dΓ − Z ∂ Ω φ w T C − µκ an e e φ dΓ+ Z ∂ Ω u w T , (cid:0) n e T21 + n e T3 s e T4 s e c T1 (cid:1) e u dΓ + Z ∂ Ω u w T , (cid:0) n e T22 + n e T3 s e T4 s e c T2 (cid:1) e u dΓ , Ω ξ T Φ − Q dΩ − Z Ω ξ T Φ − b T ε ( u ) dΩ + Z Ω ξ T , c T5 E ( φ ) dΩ + Z Ω ξ T , c T6 E ( φ ) dΩ+ Z Γ h ∪ ∂ Ω u n ξ T Φ − (cid:16) b T n e T1 − b Q T n e T (cid:17)o [[ u ]] dΓ − Z Γ h (cid:2)(cid:2) ξ T (cid:3)(cid:3) n e T6 { E ( φ ) } dΓ+ Z Γ h ∪ ∂ Ω φ (cid:8) ξ T , (cid:9) n e T51 [[ φ ]] dΓ + Z Γ h ∪ ∂ Ω φ (cid:8) ξ T , (cid:9) n e T52 [[ φ ]] dΓ − Z ∂ Ω ω ξ T , n e T6 s e s e T c T3 φ dΓ − Z ∂ Ω ω ξ T , n e T6 s e s e T c T4 φ dΓ − Z ∂ Ω Z ξ T n e T6 n e n e T E ( φ ) dΓ + Z ∂ Ω Z ξ T n e T6 n e λ Z dΓ+ X C h (cid:2)(cid:2) ξ T (cid:3)(cid:3) n e T6 s e { φ } (76) = − Z ∂ Ω P ξ T n e T6 n e e P dΓ + Z ∂ Ω u ξ T Φ − (cid:16) b T n e T1 − b Q T n e T (cid:17) e u dΓ+ Z ∂ Ω φ ξ T , (cid:0) n e T51 + n e T6 s e s e T c T3 (cid:1) e φ dΓ + Z ∂ Ω φ ξ T , (cid:0) n e T52 + n e T6 s e s e T c T4 (cid:1) e φ dΓ , Z ∂ Ω R ˆ λ T R n e n e µ dΓ = Z ∂ Ω R ˆ λ T R e R dΓ , (77) Z ∂ Ω Z ˆ λ Z n e T n e Q dΓ = Z ∂ Ω Z ˆ λ Z e Z dΓ . (78)Note that when C µκ = or Φ = , C − µκ or Φ − can be considered as an infinite large value, andthe above equations can be satisfied by letting the corresponding coefficients of C − µκ or Φ − to bezero. For example, when C µκ = , the linear coefficient of C − µκ in Eqn. 75 equals zero: − Z Ω w T µ dΩ − Z Ω w T aE ( φ ) dΩ − Z Γ h ∪ ∂ Ω φ (cid:8) w T (cid:9) an e [[ φ ]] dΓ = − Z ∂ Ω φ w T an e e φ dΓ , (79)The terms containing C − µκ in Eqn. 74 are balanced as long as Eqn. 79 is satisfied. Therefore, theycan be eliminated. Note that when C µκ = , the strain gradient κ can not be achieved directlyfrom the constitutive relation. Therefore, a post-processing algorithm for κ is required based on the36onstant terms in Eqn.75: − Z Ω w T κ dΩ − Z Ω w T , c ˆ ε ( u ) dΩ − Z Ω w T , c ˆ ε ( u ) dΩ + Z Γ h (cid:2)(cid:2) w T (cid:3)(cid:3) n e T3 { ˆ ε ( u ) } dΓ+ Z Γ h ∪ ∂ Ω u (cid:8) w T , n e T21 + w T , n e T22 (cid:9) [[ u ]] dΓ + Z ∂ Ω R w T n e T3 n e T4 n e ˆ ε ( u ) dΓ − Z ∂ Ω Q w T , n e T3 s e T4 s e c T1 u dΓ − Z ∂ Ω Q w T , n e T3 s e T4 s e c T2 u dΓ + X C h (cid:2)(cid:2) w T (cid:3)(cid:3) n e T3 s e T4 { u } + Z ∂ Ω R w T n e T3 n e T4 λ R dΓ = − Z ∂ Ω d w T n e T3 n e T4 e d dΓ+ Z ∂ Ω u w T , (cid:0) n e T21 + n e T3 s e T4 s e c T1 (cid:1) e u dΓ + Z ∂ Ω u w T , (cid:0) n e T22 + n e T3 s e T4 s e c T2 (cid:1) e u dΓ . (80)Combining Eqn. 73 – 78, we get the algebraic equations in matrix form in terms of the globalnodal vectors ( u , φ , µ , Q , λ R and λ Z ). K Muu K Muφ K Muµ K MuQ (cid:0) K Muφ (cid:1) T K Mφφ K Mφµ K MφQ (cid:0) K Muµ (cid:1) T (cid:0) K Mφµ (cid:1) T K Mµµ
MµR (cid:0) K MuQ (cid:1) T (cid:0) K MφQ (cid:1) T MQQ
MQZ (cid:0) K MµR (cid:1) T (cid:0) K MQZ (cid:1) T u φµ Q λ R λ Z = f Mu f Mφ f Mµ f MQ f MR f MZ . (81) In the reduced theory, the high-order electric displacement Q vanishes. The displacement u ,electric potential φ , high-order stress µ are still employed as independent variables at each pointand interpolated by linear trial functions. Thus the basic mixed FPM formula for the reduced flex-37electric theory can be written as: Z Ω ε T ( v ) C σε ε ( u ) dΩ − Z Ω ε T ( v ) eE ( φ ) dΩ − Z Ω ˆ ε T ( v ) c T9 µ , dΩ − Z Ω ˆ ε T ( v ) c T10 µ , dΩ − Z Γ h ∪ ∂ Ω u (cid:2)(cid:2) v T (cid:3)(cid:3) n e C σε { ε ( u ) } dΓ − Z Γ h ∪ ∂ Ω u (cid:8) ε T ( v ) (cid:9) C σε n e T1 [[ u ]] dΓ+ Z Γ h ∪ ∂ Ω u (cid:2)(cid:2) v T (cid:3)(cid:3) n e e { E ( φ ) } dΓ + Z Γ h ∪ ∂ Ω u (cid:2)(cid:2) v T (cid:3)(cid:3) { n e µ , + n e µ , } dΓ − Z Γ h ∪ ∂ Ω φ (cid:8) ε T ( v ) (cid:9) en e [[ φ ]] dΓ + Z Γ h (cid:8) ˆ ε T ( v ) (cid:9) n e [[ µ ]] dΓ − Z ∂ Ω Q v T (cid:0) c s e T4 s e µ , + c s e T4 s e µ , (cid:1) dΓ + Z ∂ Ω R ˆ ε T ( v ) n e T4 n e n e µ dΓ+ Z ∂ Ω u η h e v T u dΓ + Z Γ h η h e (cid:2)(cid:2) v T (cid:3)(cid:3) [[ u ]] dΓ + X C h (cid:8) v T (cid:9) s e n e [[ µ ]]= Z Ω v T b dΩ + Z ∂ Ω Q v T e Q dΓ − Z ∂ Ω φ ε T ( v ) en e e φ dΓ + Z ∂ Ω R ˆ ε T ( v ) n e T4 e R dΓ − Z ∂ Ω u ε T ( v ) C σε n e T1 e u dΓ + Z ∂ Ω u η h e v T e u dΓ , (82) − Z Ω E T ( τ ) (cid:0) Λ + a T C − µκ a (cid:1) E ( φ ) dΩ − Z Ω E T ( τ ) e T ε ( u ) dΩ − Z Ω E T ( τ ) a T C − µκ µ dΩ − Z Γ h ∪ ∂ Ω φ [[ τ ]] n e T (cid:0) Λ + a T C − µκ a (cid:1) { E ( φ ) } dΓ − Z Γ h ∪ ∂ Ω φ [[ τ ]] n e T e T { ε ( u ) } dΓ − Z Γ h ∪ ∂ Ω φ (cid:8) E T ( τ ) (cid:9) (cid:0) Λ + a T C − µκ a (cid:1) n e [[ φ ]] dΓ − Z Γ h ∪ ∂ Ω φ [[ τ ]] n e T a T C − µκ { µ } dΓ+ Z Γ h ∪ ∂ Ω u (cid:8) E T ( τ ) (cid:9) e T n e T1 [[ u ]] dΓ + Z ∂ Ω φ η h e τ φ dΓ + Z Γ h η h e [[ τ ]] [[ φ ]] dΓ (83) = − Z Ω τ q dΩ − Z ∂ Ω ω τ e ω dΓ − Z ∂ Ω φ E T ( τ ) (cid:0) Λ + a T C − µκ a (cid:1) n e e φ dΓ+ Z ∂ Ω u E T ( τ ) e T n e T1 e u dΓ + Z ∂ Ω φ η h e τ e φ dΓ . Z Ω w T C − µκ µ dΩ − Z Ω w T C − µκ aE ( φ ) dΩ − Z Ω w T , c ˆ ε ( u ) dΩ − Z Ω w T , c ˆ ε ( u ) dΩ+ Z Γ h ∪ ∂ Ω u (cid:8) w T , n e T21 + w T , n e T22 (cid:9) [[ u ]] dΓ − Z Γ h ∪ ∂ Ω φ (cid:8) w T (cid:9) C − µκ an e [[ φ ]] dΓ+ Z Γ h (cid:2)(cid:2) w T (cid:3)(cid:3) n e T3 { ˆ ε ( u ) } dΓ + Z ∂ Ω R w T n e T3 n e T4 n e ˆ ε ( u ) dΓ + Z ∂ Ω R w T n e T3 n e T4 λ R dΓ − Z ∂ Ω Q w T , n e T3 s e T4 s e c T1 u dΓ − Z ∂ Ω Q w T , n e T3 s e T4 s e c T2 u dΓ + X C h (cid:2)(cid:2) w T (cid:3)(cid:3) n e T3 s e T4 { u } (84) = − Z ∂ Ω d w T n e T3 n e T4 e d dΓ − Z ∂ Ω φ w T C − µκ an e e φ dΓ+ Z ∂ Ω u w T , (cid:0) n e T21 + n e T3 s e T4 s e c T1 (cid:1) e u dΓ + Z ∂ Ω u w T , (cid:0) n e T22 + n e T3 s e T4 s e c T2 (cid:1) e u dΓ , Z ∂ Ω R ˆ λ T R n e n e µ dΓ = Z ∂ Ω R ˆ λ T R e R dΓ . (85)Or in matrices form: K Muu K Muφ K Muµ (cid:0) K Muφ (cid:1) T K Mφφ K Mφµ (cid:0) K Muµ (cid:1) T (cid:0) K Mφµ (cid:1) T K Mµµ K MµR (cid:0) K MµR (cid:1) T u φµλ R = f Mu f Mφ f Mµ f MR . (86) In Eqn. 81, there are 13 independent variables at each point (3 primal variables which are thesame as in the primal method, and 10 high-order variables) and n R + n Z Lagrange multipliers.However, the number of unknown variables can be reduced by transforming the high-order variablesback to the primal displacement and electric potential at each Fragile Point at the local level. First,we rearrange Eqn. 81 as: K Mµµ K MµR (cid:0) K MµR (cid:1) T µλ R = f Mµ − (cid:0) K Muµ (cid:1) T u − (cid:0) K Mφµ (cid:1) T φ f MR , (87) K MQQ K MQZ (cid:0) K MQZ (cid:1) T Q λ Z = f MQ − (cid:0) K MuQ (cid:1) T u − (cid:0) K MφQ (cid:1) T φ f MZ . (88)Thus we can obtain a linear transformation relation between the two related nodal variable sets : µ = − G µ (cid:0) K Muµ (cid:1) T u − G µ (cid:0) K Mφµ (cid:1) T φ + g µ , (89)39 = − G Q (cid:0) K MuQ (cid:1) T u − G Q (cid:0) K MφQ (cid:1) T φ + g Q , (90)where G µ = (cid:0) K Mµµ (cid:1) − − (cid:0) K Mµµ (cid:1) − K MµR (cid:16)(cid:0) K MµR (cid:1) T (cid:0) K Mµµ (cid:1) − K MµR (cid:17) − (cid:0) K MµR (cid:1) T (cid:0) K Mµµ (cid:1) − , G Q = (cid:0) K MQQ (cid:1) − − (cid:0) K MQQ (cid:1) − K MQZ (cid:16)(cid:0) K MQZ (cid:1) T (cid:0) K MQQ (cid:1) − K MQZ (cid:17) − (cid:0) K MQZ (cid:1) T (cid:0) K MQQ (cid:1) − , g µ = G µ f Mµ + (cid:0) K Mµµ (cid:1) − K MµR (cid:16)(cid:0) K MµR (cid:1) T (cid:0) K Mµµ (cid:1) − K MµR (cid:17) − f MR , g Q = G Q f MQ + (cid:0) K MQQ (cid:1) − K MQZ (cid:16)(cid:0) K MQZ (cid:1) T (cid:0) K MQQ (cid:1) − K MQZ (cid:17) − f MZ . According to Eqn. 75 and 76, when the Fragile Points are placed at the centroid of each subdo-main, and only one integration point is used, (cid:0) K Mµµ (cid:1) − and (cid:0) K MQQ (cid:1) − are known banded matrices: (cid:0) K Mµµ (cid:1) − = − diag (cid:18)(cid:20) A , A , · · · , A N (cid:21)(cid:19) ⊗ C µκ , (cid:0) K MQQ (cid:1) − = diag (cid:18)(cid:20) A , A , · · · , A N (cid:21)(cid:19) ⊗ Φ , (91)where A i denotes the area of each subdomain E i , and N is the number of all the Fragile Points.Therefore, the transforming matrices G µ and G Q can be easily achieved by inverting a n R × n R matrix h(cid:0) K MµR (cid:1) T (cid:0) K Mµµ (cid:1) − K MµR i and a n Z × n Z matrix h(cid:0) K MQZ (cid:1) T (cid:0) K MQQ (cid:1) − K MQZ i respectively.We can prove that G µ and G Q are also banded. All the Fragile Points can be divided into twogroups: for P i ∈ b E in , the nodal displacement or electric potential value at P i has no influence ontrial functions on the external boundary; while at P i ∈ b E ex , a non-zero nodal value will result ina non-zero trial function on ∂ Ω , i.e., b E ex = S ∂E i ∪ ∂ Ω = ∅ b E i , where b E i is the set of all the supportingpoints of P i . Thus, we have: µ i = G iµu u iE + G iµφ φ iE , for P i ∈ b E in G iµu u exE + G iµφ φ exE + g iµ , for P i ∈ b E ex , Q i = G iQu u iE + G iQφ φ iE , for P i ∈ b E in G iQu u exE + G iQφ φ exE + g iQ , for P i ∈ b E ex , (92)where µ i and Q i are the nodal values of µ and Q at P i . u iE and φ iE are the nodal displacementand electric potential vectors for all P j ∈ b E i , where b E i = S P k ∈ b E i b E k is the “generalizedâĂŹâĂŹ setof supporting points of P i . u exE and φ exE are the nodal displacement and electric potential vectors40or all P j ∈ b E exE , where b E exE = S P k ∈ b E ex b E k . G iµu , G iµφ , G iQu , G iQφ , g iµ and g iQ are linear algebraicmatrices (vectors).Thus, all the high-order variables are eliminated at the point level. Finally, the mixed FPMformula with only nodal displacement and electric potential is achieved: K Muu K Muφ (cid:16) K Muφ (cid:17) T K Mφφ u φ = f Mu f Mφ , or K M x = f M . (93)In full flexoelectric theory: K Muu = K Muu − K Muµ G µ (cid:0) K Muµ (cid:1) T − K MuQ G Q (cid:0) K MuQ (cid:1) T , K Muφ = K Muφ − K Muµ G µ (cid:0) K Mφµ (cid:1) T − K MuQ G Q (cid:0) K MφQ (cid:1) T , K Mφφ = K Mφφ − K Mφµ G µ (cid:0) K Mφµ (cid:1) T − K MφQ G Q (cid:0) K MφQ (cid:1) T , f Mu = f Mu − K Muµ g µ − K MuQ g Q , f Mφ = f Mφ − K Mφµ g µ − K MφQ g Q , While in the reduced theory, the last terms relating to Q in the above equations are absent. Thesame as in the primal FPM, when the full flexoelectric theory is employed, Eqn. 93 is a nonlinearequation system and can be solved by Newton-Raphson method. The stiffness matrix K M is sparseand symmetric for both the theories. Also, there are only three explicit unknown variables at eachFragile Point.
5. Simulations of Crack Initiation & Propagation in FPM
As a benefit of the discontinuous trial and test functions, it is much simpler in the FPM tosimulate crack initiation and propagation, as compared with other numerical methods. In the presentwork, all the cracks are considered as mechanically traction-free. Yet in the electrical field, twodifferent crack-face boundary conditions are considered. The first one is known as an electricallyimpermeable boundary condition, in which: e ω + = e ω − = 0 , ∆ φ = φ + − φ − = 0 . (94)This condition implies that the gap between the upper ( + ) and lower ( − ) crack-faces are filled with amedium with zero electric permittivity. Alternatively, an electrically permeable crack-face boundary41ondition is defined as: φ + = φ − , e ω + = e ω − , (95)that is, an infinite electric permittivity of the medium between the upper and lower crack-faces isassumed.When an electrically impermeable crack emerges between two adjacent subdomains, we convertthe corresponding internal boundary to two external boundaries with zero mechanical and electricaltractions ( e ∈ ∂ Ω Q ∩ ∂ Ω R ∩ ∂ Ω ω ∩ ∂ Ω Z ). And we also cut off the interaction between the two Pointsin the adjacent subdomains, i.e., the Points will be removed from the neighboring point set of eachother. Therefore, the support of the point on one side of the crack no longer contain the adjacentpoints and its neighboring points on the other side of the crack. For example, in Fig. 4, as a resultof the crack (presented by a red polyline), three points (grey) are removed from the support of thereference point (red). Figure 4: Support of a point by the sides of cracks.
As for the electrically permeable crack-face boundaries, two different interpolations are em-ployed in the mechanical and electrical field. In the mechanical field, the internal boundaries on thecracks are converted into traction-free boundaries and interactions between the points on either side42f the crack are cut off. Whereas in the electrical interpolation, these boundaries remain as internalboundaries with continuous electrical field. As a result, the support of points and definition of Γ h may be different in the mechanical and electrical interpolations.Therefore, for both the crack-face boundaries conditions, only slight adjustments are required forthe stiffness matrix and/or load vectors when a crack occurs. That is, for the two points adjacent tothe crack and their closest neighbors, the point stiffness matrices will be regenerated, and the termsrelated to the numerical fluxes on the cracked boundaries should be deleted. The total numberof DoF remains the same. And in the linear reduced theory, the load vector does not need to beadjusted. This is much more convenient than remeshing or deleting elements as in the FEM whensimulating the crack propagation.Various crack initiation and propagation criteria can be applied and incorporated within theFPM. The choice of these criteria should be based on the characteristics of the realistic problemand material. In this paper, for example, the classic Maximum Hoop Stress criterion and an Inter-Subdomain-Boundary Bonding-Energy-Rate(BER)-Based criterion are employed.In crack initiation analysis, we consider a quasi-static loading process. The entire load historycan be divided into several steps. And in each step, for simplicity, the most “dangerous” internalboundaries will be cracked. That is, for instance, based on the Maximum Hoop Stress criterion, theinternal boundary with the largest normal traction will be cracked. Alternatively, here we introducean energy-based criterion initiated by the J-integral. Based on the electric Gibbs free energy G , theJ-integral in domain Ω int is defined as [55]: J = Z ∂ Ω int Gn dΓ − Z ∂ Ω int (cid:16) ∂ u ∂ ˆ x (cid:17) T e Q + (cid:16) ∂ ( n e ˆ ε ( u )) ∂ ˆ x (cid:17) T e R − ∂φ∂ ˆ x e ω − (cid:18) ∂ ( n e T E ( φ ) ) ∂ ˆ x (cid:19) T e Z dΓ . (96)The previous equation is written in a local ˆ x - ˆ x coordinate system, in which ˆ x is aligned with thecrack and ˆ x is aligned with the normal vector. e Q , e R , e ω and e Z are defined in Eqn. 10 – 13, yet theyare unknown at here. When Ω int is divided into several subdomains Ω intI , we can rewrite Eqn. 9643s: J = X I Z Ω intI (cid:20) ∂ ˆ ε ( u ) ∂ ˆ x ∂G∂ ˆ ε + ∂ κ ( u ) ∂ ˆ x ∂G∂ κ + ∂ E ( φ ) ∂ ˆ x ∂G∂ E + ∂ V ( φ ) ∂ ˆ x ∂G∂ V (cid:21) dΩ − Z ∂ Ω int "(cid:18) ∂ u ∂ ˆ x (cid:19) T e Q + ∂ ˆ ε T ( u ) ∂ ˆ x n e T4 e R − ∂φ∂ ˆ x e ω − ∂ E T ( φ ) ∂ ˆ x n e e Z dΓ − X Z Γ inth [[ Gn e ]] dΓ= X I Z Ω intI (cid:20) ∂ ε ( u ) ∂ ˆ x σ + ∂ ˆ ε ( u ) ∂ ˆ x σ ES + ∂ κ ( u ) ∂ ˆ x µ − ∂ E ( φ ) ∂ ˆ x D + ∂ V ( φ ) ∂ ˆ x Q (cid:21) dΩ − Z ∂ Ω int "(cid:18) ∂ u ∂ ˆ x (cid:19) T e Q + ∂ ˆ ε T ( u ) ∂ ˆ x n e T4 e R − ∂φ∂ ˆ x e ω − ∂ E T ( φ ) ∂ ˆ x n e e Z dΓ − X Z Γ inth [[ Gn e ]] dΓ (97)where Γ inth is the set of internal boundaries within Ω int . Since ˆ x is aligned with the crack, n e = 0 thus the last term in Eqn. 97 can be eliminated.The FPM formulations Eqn. 56 and 57 should be satisfied in the present domain Ω int with anyarbitrary test functions v and τ . Therefore, here we let v = ∂ u ∂ ˆ x and τ = ∂φ∂ ˆ x and comparing theFPM equations with Eqn 97. The essential boundaries are omitted here. Finally, we can achieve: J ≈ X e ∈ Γ inth BER, (98)44here
BER = − Z e η h e (cid:20)(cid:20) ∂ u T ∂ ˆ x (cid:21)(cid:21) [[ u ]] dΓ − Z e η h e (cid:20)(cid:20) ∂ ˆ ε T ( u ) ∂ ˆ x (cid:21)(cid:21) n e T4 n e [[ ˆ ε ( u )]] dΓ − Z e η h e (cid:20)(cid:20) ∂φ∂ ˆ x (cid:21)(cid:21) [[ φ ]] dΓ − Z e η h e (cid:20)(cid:20) ∂ E T ( φ ) ∂ ˆ x n e (cid:21)(cid:21) n e T4 n e (cid:2)(cid:2) n e T E ( φ ) (cid:3)(cid:3) dΓ+ Z e (cid:18)(cid:20)(cid:20) ∂ u T ∂ ˆ x (cid:21)(cid:21) n e C σε { ε ( u ) } + (cid:26) ∂ ε T ( u ) ∂ ˆ x (cid:27) C σε n e T1 [[ u ]] (cid:19) dΓ+ Z e (cid:18)(cid:20)(cid:20) ∂ ˆ ε T ( u ) ∂ ˆ x (cid:21)(cid:21) n e C µκ { κ ( u ) } + (cid:26) ∂ κ T ( u ) ∂ ˆ x (cid:27) C µκ n e T3 [[ ˆ ε ( u )]] (cid:19) dΓ+ Z e (cid:18)(cid:20)(cid:20) ∂φ∂ ˆ x (cid:21)(cid:21) n e T Λ { E ( φ ) } + (cid:26) ∂ E T ( φ ) ∂ ˆ x (cid:27) Λn e [[ φ ]] (cid:19) dΓ − Z e (cid:18)(cid:20)(cid:20) ∂ E T ( φ ) ∂ ˆ x (cid:21)(cid:21) n e Φ { V ( φ ) } + (cid:26) ∂ V T ( φ ) ∂ ˆ x (cid:27) Φn e T6 [[ E ( φ )]] (cid:19) dΓ − Z e (cid:20)(cid:20) ∂ u T ∂ ˆ x (cid:21)(cid:21) n e C µκ { κ , ( u ) } dΓ − Z e (cid:20)(cid:20) ∂ u T ∂ ˆ x (cid:21)(cid:21) n e C µκ { κ , ( u } ) dΓ − Z e (cid:20)(cid:20) ∂ u T ∂ ˆ x (cid:21)(cid:21) n(cid:16) n e b − n e ac − n e ac − n e e Q ( u , φ ) (cid:17) V ( φ ) o dΓ − Z e (cid:26) ∂ V T ( φ ) ∂ ˆ x (cid:16) b T n e T1 − c T5 a T n e T21 − c T6 a T n e T22 − e Q T ( u , φ ) n e T4 (cid:17)(cid:27) [[ u ]] dΓ − Z e (cid:20)(cid:20) ∂ u T ∂ ˆ x (cid:21)(cid:21) n(cid:16) n e e − n e e D ( u , φ ) (cid:17) E ( φ ) o dΓ + Z e (cid:26) ∂ ε T ( u ) ∂ ˆ x (cid:27) en e [[ φ ]] dΓ − Z e (cid:26) ∂ E T ( φ ) ∂ ˆ x (cid:16) e T n e T1 − e D T ( u , φ ) n e T4 (cid:17)(cid:27) [[ u ]] dΓ + Z e (cid:20)(cid:20) ∂φ∂ ˆ x (cid:21)(cid:21) n e T e T { ε ( u ) } dΓ+ Z e (cid:26) ∂ κ T ( u ) ∂ ˆ x (cid:27) (cid:0) an e − c T7 bn e T51 − c T8 bn e T52 (cid:1) [[ φ ]] dΓ+ Z e (cid:20)(cid:20) ∂φ∂ ˆ x (cid:21)(cid:21) (cid:0) n e T a T − n e b T c − n e b T c (cid:1) { κ ( u ) } dΓ − Z e (cid:20)(cid:20) ∂ ˆ ε T ( u ) ∂ ˆ x (cid:21)(cid:21) n e a { E ( φ ) } dΓ − Z e (cid:26) ∂ ε T ( u ) ∂ ˆ x (cid:27) bn e T6 [[ E ( φ )]] dΓ − Z e (cid:26) ∂ E T ( φ ) ∂ ˆ x (cid:27) a T n e T3 [[ ˆ ε ( u )]] dΓ − Z e (cid:20)(cid:20) ∂ E T ( φ ) ∂ ˆ x (cid:21)(cid:21) n e b T { ε ( u ) } dΓ − Z e (cid:18)(cid:20)(cid:20) ∂φ∂ ˆ x (cid:21)(cid:21) n e Φ { V , ( φ ) } (cid:19) dΓ − Z e (cid:18)(cid:20)(cid:20) ∂φ∂ ˆ x (cid:21)(cid:21) n e Φ { V , ( φ ) } (cid:19) dΓ . (99)Thus, the J-integral can be approximated by a sum of integrals over all the internal boundaries.Here we postulate the component on each internal boundary as an estimate of the bonding energy rate( BER ). In the current work, the
BER is employed as an energy-based criterion of crack initiationand development. An internal boundary segment will be cracked if its
BER exceeds a prescribed45ritical value. Note that in the FPM, the
BER on each internal boundary is generated following thesame process of the internal boundary stiffness matrix K h , except that the test functions v and τ arereplaced by ∂ u ∂ ˆ x and ∂φ∂ ˆ x , where u and φ are solutions of the FPM analysis.
6. Conclusion
A meshless Fragile Points Method (FPM) based on Galerkin weak form is developed to analyzeflexoelectric effects in dielectric solids. Both primal and mixed formulations are presented. Local,simple, polynomial and discontinuous trial and test functions, which take the form of Taylor expan-sions at each internal Fragile Point, are employed. Local RBF-DQ method is applied to approximatethe higher derivatives. Numerical Flux Corrections are ultilized to ensure the continuity condition.In the mixed FPM, all the additional high order independent variables are eliminated locally at eachPoint, and only 3 DoFs are retained explicitly in the final formula. Both full and reduced flexoelec-tric theories are taken into consideration. The full theory involving the electrostatic stress leads anonlinear algebraic system while the formulations for the reduced theory are linear. The algorithmfor simulating crack initiation and propagation in the FPM is also presented, using the stress-basedcriterion as well as an Inter-Subdomain-Boundary Bonding-Energy-Rate(BER)-Based criterion forcrack initiation and development.To conclude, the primal and mixed FPM proposed in this work have the following features andadvantages as compared to other computational methods for flexoelectricity published in prior lit-erature:• The FPM is a meshless or element-free method with point-based trial and test functions thatare piece-wise continuous.• The integration of the Galerkin weak form is very simple by using Guassian quadrature withonly one integration point being sufficient most of the time.• The trial function has the Kronecker-delta property; the essential boundary conditions can beimposed directly or weakly using numerical fluxes.• Crack and rupture initiation and propagation can be easily simulated by FPM and do notinvolve remeshing or trial function enhancement.46 Crack and rupture criteria are based on simple and commonly used continuum physics.• The stiffness matrix is sparse and symmetric.• With mixed FPM, the first and higher strain gradient problems can be handled with ease. Onlythe lowest-order primitive variables are retained as DoFs at each node.• There is no incompressibility locking or shear locking in FPM.• The FPM does not suffer from mesh distortion.In this first part of the present two-paper series, the theoretical formulation and implementationof the proposed methodology are given in detail. Numerical examples and validation are reportedin Part II of the present study.
References [1] B. Wang, Y. Gu, S. Zhang, and L. Q. Chen. Flexoelectricity in solids: Progress, challenges,and perspectives.
Progress in Materials Science , 106(September 2018), 2019.[2] X. Zhuang, B. H. Nguyen, S. S. Nanthakumar, T. Q. Tran, N. Alajlan, and T. Rabczuk. Com-putational modeling of flexoelectricity-A review.
Energies , 16(3):1–30, 2020.[3] R. Maranganti, N. D. Sharma, and P. Sharma. Electromechanical coupling in nonpiezoelectricmaterials due to nanoscale nonlocal size effects: Green’s function solutions and embeddedinclusions.
Physical Review B - Condensed Matter and Materials Physics , 74(1):1–14, 2006.[4] W. Huang, X. Yan, S. R. Kwon, S. Zhang, F. G. Yuan, and X. Jiang. Flexoelectric straingradient detection using Ba0.64Sr 0.36TiO3 for sensing.
Applied Physics Letters , 101(25):0–4, 2012.[5] S.-B. Choi and G.-W. Kim. Measurement of flexoelectric response in polyvinylidene fluoridefilms for piezoelectric vibration energy harvesters.
Journal of Physics D: Applied Physics , 50(7):75502, 2017. 476] C. Liu, J. Wang, G. Xu, M. Kamlah, and T. Y. Zhang. An isogeometric approach to flexoelectriceffect in ferroelectric materials.
International Journal of Solids and Structures , 162:198–210,2019.[7] S. Shen and S. Hu. A theory of flexoelectricity with surface effect for elastic dielectrics.
Journal of the Mechanics and Physics of Solids , 58(5):665–677, 2010.[8] S. Hu and S. Shen. Variational principles and governing equations in nano-dielectrics with theflexoelectric effect.
Science China: Physics, Mechanics and Astronomy , 53(8):1497–1504,2010.[9] M. S. Majdoub, P. Sharma, and T. Cagin. Enhanced size-dependent piezoelectricity and elas-ticity in nanostructures due to the flexoelectric effect.
Phys. Rev. B , 77(12):125424, 2008.[10] S. Mao and P. K. Purohit. Insights into flexoelectric solids from strain-gradient elasticity.
Journal of Applied Mechanics, Transactions ASME , 81(8):1–10, 2014.[11] J. Yvonnet and L. P. Liu. A numerical framework for modeling flexoelectricity and Maxwellstress in soft dielectrics at finite strains.
Computer Methods in Applied Mechanics and Engi-neering , 313:450–482, 2017.[12] S. Mao, P. K. Purohit, and N. Aravas. Mixed finite-element formulations in piezoelectricity andflexoelectricity.
Proceedings of the Royal Society A: Mathematical, Physical and EngineeringSciences , 472(2190), 2016.[13] J. Sladek, V. Sladek, M. W¨unsche, and C. Zhang. Effects of electric field and strain gradientson cracks in piezoelectric solids.
European Journal of Mechanics, A/Solids , 71(March):187–198, 2018.[14] F. Deng, Q. Deng, and S. Shen. A Three-Dimensional Mixed Finite Element for Flexoelec-tricity.
Journal of Applied Mechanics, Transactions ASME , 85(3), 2018.[15] B. He, B. Javvaji, and X. Zhuang. Characterizing flexoelectricity in composite material usingthe element-free Galerkin method.
Energies , 12(2), 2019.4816] J. Sladek, V. Sladek, and P. H. Wen. The meshless analysis of scale-dependent problems forcoupled fields.
Materials , 13(11), 2020.[17] A. Beheshti. Finite element analysis of plane strain solids in strain-gradient elasticity.
ActaMechanica , 228(10):3543–3559, 2017.[18] J. Sladek, V. Sladek, M. W¨unsche, and C. L. Tan. Flexoelectric effect for cracks in piezoelectricsolids.
Key Engineering Materials , 774 KEM:90–95, 2018.[19] J. Y. Shu, W. E. King, and N. A. Fleck. Finite elements for materials with strain gradienteffects.
International Journal for Numerical Methods in Engineering , 44(3):373–391, 1999.[20] E. Amanatidou and N. Aravas. Mixed finite element formulations of strain-gradient elastic-ity problems.
Computer Methods in Applied Mechanics and Engineering , 191(15-16):1723–1751, 2002.[21] P. L. Bishay, J. Sladek, V. Sladek, and S. N. Atluri. Analysis of functionally graded magneto-electro-elastic composites using hybrid/mixed finite elements and node-wise material proper-ties.
Computers, Materials and Continua , 29(3):213–261, 2012.[22] F. Deng, Q. Deng, W. Yu, and S. Shen. Mixed finite elements for flexoelectric solids.
Journalof Applied Mechanics, Transactions ASME , 84(8), 2017.[23] Y. Mao, S. Ai, X. Xiang, and C. Chen. Theory for dielectrics considering the direct andconverse flexoelectric effects and its finite element implementation.
Applied MathematicalModelling , 40(15-16):7115–7137, 2016.[24] S. S. Nanthakumar, X. Zhuang, H. S. Park, and T. Rabczuk. Topology optimization of flexo-electric structures.
Journal of the Mechanics and Physics of Solids , 105:217–234, 2017.[25] A. Abdollahi, C. Peco, D. Mill´an, M. Arroyo, and I. Arias. Computational evaluation of theflexoelectric effect in dielectric solids.
Journal of Applied Physics , 116(9), 2014.[26] M. C. Ray. Mesh-free models for static analysis of smart laminated composite beams.
Inter-national Journal for Numerical Methods in Engineering , 109:1804–1820, 2017.4927] M. C. Ray. Mesh free model of nanobeam integrated with a flexoelectric actuator layer.
Com-posite Structures , 159:63–71, 2017.[28] R. Basutkar, S. Sidhardh, and M. C. Ray. Static analysis of flexoelectric nanobeams incorpo-rating surface effects using element free Galerkin method.
European Journal of Mechanics,A/Solids , 76(December 2018):13–24, 2019.[29] S. Sidhardh and M. C. Ray. Element-free Galerkin model of nano-beams considering straingradient elasticity.
Acta Mechanica , 229(7):2765–2786, 2018.[30] Z. Tang, S. Shen, and S. N. Atluri. Analysis of materials with strain-gradient effects: A mesh-less local Petrov-Galerkin (MLPG) approach, with nodal displacements only.
CMES - Com-puter Modeling in Engineering and Sciences , 4(1):177–196, 2003.[31] J. Sladek, V. Sladek, P. Stanak, P. H. Wen, and S. N. Atluri. Laminated elastic plates withpiezoelectric sensors and actuators.
CMES - Computer Modeling in Engineering and Sciences ,85(6):543–572, 2012.[32] S. N. Atluri, Z. D. Han, and A. M. Rajendran. A new implementation of the meshless fi-nite volume method, through the MLPG ”mixed” approach.
CMES - Computer Modeling inEngineering and Sciences , 6(6):491–513, 2004.[33] S. N. Atluri and S. Shen. Simulation of a 4th order ODE: Illustration of various primal &mixed MLPG methods.
Laser Physics , 15(3):241–268, 2005.[34] T. Jarak, B. Jaluˇsi´c, and J. Sori´c. Mixed meshless local Petrov-Galerkin methods for solvinglinear fourth-order differential equations.
Transactions of Famena , 44(1):1–12, 2020.[35] X. Zhuang, S. S. Nanthakumar, and T. Rabczuk. A meshfree formulation for large defor-mation analysis of flexoelectric structures accounting for the surface effects. arXiv preprintarXiv:1911.06553 , 2019.[36] T. Zhu and S. N. Atluri. A modified collocation method and a penalty formulation for enforc-ing the essential boundary conditions in the element free Galerkin method.
ComputationalMechanics , 21(3):211–222, 1998. 5037] M. Hillman and J. Chen. An accelerated, convergent, and stable nodal integration in Galerkinmeshfree methods for linear and nonlinear mechanics.
International Journal for NumericalMethods in Engineering , 107(7):603–630, 2016.[38] Q. Duan, X. Li, H. Zhang, and T. Belytschko. SecondâĂŘorder accurate derivatives andintegration schemes for meshfree methods.
International Journal for Numerical Methods inEngineering , 92(4):399–424, 2012.[39] L. Dong, T. Yang, K. Wang, and S. N. Atluri. A new Fragile Points Method (FPM) in computa-tional mechanics, based on the concepts of Point Stiffnesses and Numerical Flux Corrections.
Engineering Analysis with Boundary Elements , 107(July):124–133, 2019.[40] Y. Guan, R. Grujicic, X. Wang, L. Dong, and S. N. Atluri. A new meshless âĂIJfragile pointsmethodâĂİ and a local variational iteration method for general transient heat conduction inanisotropic nonhomogeneous media. Part I: Theory and implementation.
Numerical HeatTransfer, Part B: Fundamentals , 78(2):71–85, aug 2020.[41] Y. Guan, R. Grujicic, X. Wang, L. Dong, and S. N. Atluri. A new meshless âĂIJfragilepoints methodâĂİ and a local variational iteration method for general transient heat conductionin anisotropic nonhomogeneous media. Part II: Validation and discussion.
Numerical HeatTransfer, Part B: Fundamentals , 78(2):86–109, aug 2020.[42] T. Yang, L. Dong, and S. N. Atluri. A Simple Galerkin Meshless Method, the Fragile PointsMethod (FPM) Using Point Stiffness Matrices, for 2D Linear Elastic Problems in ComplexDomains with Crack and Rupture Propagation. arXiv , pages arXiv–1909, 2019.[43] H. V. Do, T. Lahmer, X. Zhuang, N. Alajlan, H. Nguyen-Xuan, and T. Rabczuk. An isogeo-metric analysis to identify the full flexoelectric complex material properties based on electricalimpedance curve.
Computers and Structures , 214:1–14, 2019.[44] V. P. Nguyen, C. Anitescu, S. P. Bordas, and T. Rabczuk. Isogeometric analysis: An overviewand computer implementation aspects.
Mathematics and Computers in Simulation , 117:89–116, 2015. 5145] J. Sladek, V. Sladek, and M. Jus. The MLPG for crack analyses in composites with flexoelec-tricity effects.
Composite Structures , 204(April):105–113, 2018.[46] A. Abdollahi and I. Arias. Phase-field modeling of crack propagation in piezoelectric andferroelectric materials with different electromechanical crack conditions.
Journal of the Me-chanics and Physics of Solids , 60(12):2100–2126, 2012.[47] A. Abdollahi, C. Peco, D. Mill´an, M. Arroyo, G. Catalan, and I. Arias. Fracture tougheningand toughness asymmetry induced by flexoelectricity.
Physical Review B - Condensed Matterand Materials Physics , 92(9), 2015.[48] P. O’Hara, J. Hollkamp, C. A. Duarte, and T. Eason. A two-scale generalized finite elementmethod for fatigue crack propagation simulations utilizing a fixed, coarse hexahedral mesh.
Computational Mechanics , 57(1):55–74, 2016.[49] T. P. Fries and T. Belytschko. The extended/generalized finite element method: An overview ofthe method and its applications.
International Journal for Numerical Methods in Engineering ,84(3):253–304, 2010.[50] B. Giovanardi, A. Scotti, and L. Formaggia. A hybrid XFEM âĂŞPhase field (Xfield) methodfor crack propagation in brittle elastic materials.
Computer Methods in Applied Mechanicsand Engineering , 320:396–420, 2017.[51] J. Hou, M. Goldstraw, S. Maan, and M. Knop. An evaluation of 3D crack growth using ZEN-CRACK. Technical report, 2001.[52] C. Shu, H. Ding, and K. S. Yeo. Local radial basis funcion-based differential quadraturemethod and its application to solve two-dimensional incompressible Navier-Stokes equations.
Computer Methods in Applied Mechanics and Engineering , 192(7-8):941–954, 2003.[53] F. Bassi and S. Rebay. A high-order accurate discontinuous finite element method for thenumerical solution of the compressible Navier-Stokes equations.
Journal of ComputationalPhysics , 131(2):267–279, 1997. 5254] I. Mozolevski, E. S¨uli, and P. R. B¨osing. hp-version a priori error analysis of interior penaltydiscontinuous Galerkin finite element approximations to the biharmonic equation.
Journal ofScientific Computing , 30(3):465–491, 2007.[55] J. Sladek, V. Sladek, P. Stanak, C. Zhang, and C. L. Tan. Fracture mechanics analysis of size-dependent piezoelectric solids.
International Journal of Solids and Structures , 113-114:1–9,2017.
Acknowledgment
Yue Guan thankfully acknowledges the financial support for her work, provided through thefunding for Professor AtluriâĂŹs Presidential Chair at TTU.53 ppendix A. Matrices in numerical implementation c = , c = c = h i , c = h i c = c ⊗ I × , c = c ⊗ I × , c = /
20 0 0 1 1 / , c = / / c = , c = n = h n n i T , s = h s s i T = ± h n − n i T n = n n n n , n = n n n n , n = n n
00 0 0 n n , n = n n n n , n = n n n n n n
00 0 0 n n , n = n n n n , s = s s s s , n = h n n i , n = h n n i ,,