Access the full text.
Sign up today, get DeepDyve free for 14 days.
A. Abdollahi, C. Peco, Daniel Millán, M. Arroyo, I. Arias (2014)
Computational evaluation of the flexoelectric effect in dielectric solidsJournal of Applied Physics, 116
J. Shu, W. King, N. Fleck (1999)
FINITE ELEMENTS FOR MATERIALS WITH STRAIN GRADIENT EFFECTSInternational Journal for Numerical Methods in Engineering, 44
Tran Thai, T. Rabczuk, X. Zhuang (2018)
A large deformation isogeometric approach for flexoelectricity and soft materialsComputer Methods in Applied Mechanics and Engineering
R. Maranganti, N. Sharma, P. Sharma (2006)
Electromechanical coupling in nonpiezoelectric materials due to nanoscale nonlocal size effects: Green's function solutions and embedded inclusionsPhysical Review B, 74
Shuwen Zhang, Kaiyuan Liu, Minglong Xu, S. Shen (2017)
A curved resonant flexoelectric actuatorApplied Physics Letters, 111
Ruijian Zhu, Zengmei Wang, He Ma, G. Yuan, Fengxia Wang, Zhenxiang Cheng, H. Kimura (2018)
Poling-free energy harvesters based on robust self-poled ferroelectric fibersNano Energy
J.-T. Liu, B. He, S. Gu, Q. He (2018)
Implementation of a physics‐based general elastic imperfect interface model in the XFEM and LSM contextInternational Journal for Numerical Methods in Engineering, 115
Xiaoning Jiang, Wenbin Huang, Shujun Zhang (2013)
Flexoelectric nano-generator: Materials, structures and devicesNano Energy, 2
Liping Liu (2013)
On energy formulations of electrostatics for continuum mediaJournal of The Mechanics and Physics of Solids, 61
R. Guo, Lei Shen, Han Wang, Z. Lim, Wenlai Lu, P. Yang, Ariando, A. Gruverman, T. Venkatesan, Yuanping Feng, Jingsheng Chen (2016)
Tailoring Self‐Polarization of BaTiO3 Thin Films by Interface Engineering and Flexoelectric EffectAdvanced Materials Interfaces, 3
S. Gu, Q. He (2011)
Interfacial discontinuity relations for coupled multifield phenomena and their application to the modeling of thin interphases as imperfect interfacesJournal of The Mechanics and Physics of Solids, 59
S. Kwon, W. Huang, S. Zhang, F. Yuan, X. Jiang (2016)
Study on a flexoelectric microphone using barium strontium titanateJournal of Micromechanics and Microengineering, 26
S. Krichen, P. Sharma (2016)
Flexoelectricity: A Perspective on an Unusual Electromechanical CouplingJournal of Applied Mechanics, 83
Longlong Shu, Xiaoyong Wei, Ting-Tian Pang, Xi Yao, Chunlei Wang (2011)
Symmetry of flexoelectric coefficients in crystalline mediumJournal of Applied Physics, 110
(1968)
Changes in curvature of a ferroelectric film due to polarization
K. Zhenbang (2007)
Some problems in electrostrictive and magnetostrictive materialsActa Mechanica Solida Sinica, 20
U. Bhaskar, N. Banerjee, A. Abdollahi, Zhe Wang, D. Schlom, G. Rijnders, G. Catalán (2016)
A flexoelectric microelectromechanical system on silicon.Nature nanotechnology, 11 3
E. Sahin, S. Dost (1988)
A strain-gradients theory of elastic dielectrics with spatial dispersionInternational Journal of Engineering Science, 26
Wenbin Huang, Xiang Yan, S. Kwon, Shujun Zhang, F. Yuan, Xiaoning Jiang (2012)
Flexoelectric strain gradient detection using Ba0.64Sr0.36TiO3 for sensingApplied Physics Letters, 101
Y. Krongauz, T. Belytschko (1998)
EFG approximation with discontinuous derivativesInternational Journal for Numerical Methods in Engineering, 41
M. Majdoub, P. Sharma, T. Çagin (2008)
Dramatic enhancement in energy harvesting for a narrow range of dimensions in piezoelectric nanostructuresPhysical Review B, 78
H. Lu, C. Bark, D. Ojos, J. Alcalá, C. Eom, G. Catalán, A. Gruverman
Supplementary Materials for Mechanical Writing of Ferroelectric Polarization Materials and Methods Supplementary Text Figs. S1 to S4 References
Q. Deng, M. Kammoun, A. Erturk, P. Sharma (2014)
Nanoscale flexoelectric energy harvestingInternational Journal of Solids and Structures, 51
A. Rey, P. Servio, E. Herrera-Valencia (2013)
Bioinspired model of mechanical energy harvesting based on flexoelectric membranes.Physical review. E, Statistical, nonlinear, and soft matter physics, 87 2
M. Majdoub, P. Sharma, T. Çagin (2009)
Erratum: Dramatic enhancement in energy harvesting for a narrow range of dimensions in piezoelectric nanostructures [Phys. Rev. B 78, 121407(R) (2008)]Physical Review B, 79
M. Majdoub, P. Sharma, T. Çagin (2009)
Erratum: Enhanced size-dependent piezoelectricity and elasticity in nanostructures due to the flexoelectric effect [Phys. Rev. B 77 , 125424 (2008)]Physical Review B, 79
Ming-Min Yang, Dong Kim, M. Alexe (2018)
Flexo-photovoltaic effectScience, 360
V. Merupo, B. Guiffard, R. Seveno, M. Tabellout, A. Kassiba (2017)
Flexoelectric response in soft polyurethane films and their use for large curvature sensingJournal of Applied Physics, 122
(1964)
Piezoelectric effect during inhomogeneous deformation and acoustic scattering of carriers in crystals
L. Piegl, W. Tiller (1995)
The NURBS Book
B. Schrefler, S. Secchi, L. Simoni (2006)
On adaptive refinement techniques in multi-field problems including cohesive fractureComputer Methods in Applied Mechanics and Engineering, 195
H. Quang, Q. He (2011)
The number and types of all possible rotational symmetries for flexoelectric tensorsProceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 467
Z. Kuang (2009)
Internal energy variational principles and governing equations in electroelastic analysisInternational Journal of Solids and Structures, 46
Vinh Nguyen, C. Anitescu, S. Bordas, T. Rabczuk (2012)
Isogeometric analysis: An overview and computer implementation aspectsArXiv, abs/1205.2129
Guirong Liu (2009)
Meshfree Methods: Moving Beyond the Finite Element Method, Second Edition
R. Borst (2017)
Fluid flow in fractured and fracturing porous media: A unified viewMechanics Research Communications, 80
N. Sukumar, Zhenyu Huang, J. Prévost, Z. Suo (2004)
Partition of unity enrichment for bimaterial interface cracksInternational Journal for Numerical Methods in Engineering, 59
B. Nguyen, X. Zhuang, T. Rabczuk (2018)
Numerical model for the characterization of Maxwell-Wagner relaxation in piezoelectric and flexoelectric composite materialComputers & Structures
H. Ghasemi, Harold Park, T. Rabczuk (2017)
A level-set based IGA formulation for topology optimization of flexoelectric materialsComputer Methods in Applied Mechanics and Engineering, 313
Haidong Lu, Shi Liu, Z. Ye, S. Yasui, H. Funakubo, A. Rappe, A. Gruverman (2017)
Asymmetry in mechanical polarization switchingApplied Physics Letters, 110
Seung-bok Choi, Gi-Woo Kim (2017)
Measurement of flexoelectric response in polyvinylidene fluoride films for piezoelectric vibration energy harvestersJournal of Physics D: Applied Physics, 50
A. Abdollahi, N. Domingo, I. Arias, G. Catalán (2019)
Converse flexoelectricity yields large piezoresponse force microscopy signals in non-piezoelectric materialsNature Communications, 10
Yucong Liu, Jiadong Chen, H. Deng, G. Hu, Da-ming Zhu, N. Dai (2016)
Anomalous thermoelectricity in strained Bi2Te3 filmsScientific Reports, 6
Liping Liu (2014)
An energy formulation of continuum magneto-electro-elasticity with applicationsJournal of The Mechanics and Physics of Solids, 63
Feng Deng, Q. Deng, Wenshan Yu, S. Shen (2017)
Mixed finite elements for flexoelectric solidsJournal of Applied Mechanics, 84
P. Zubko, G. Catalán, A. Tagantsev (2013)
Flexoelectric Effect in SolidsAnnual Review of Materials Research, 43
(1963)
Long wavelength oscillations of diamond-type crystals including long range forces
N. Sukumar, D. Chopp, N. Moës, T. Belytschko (2001)
MODELING HOLES AND INCLUSIONS BY LEVEL SETS IN THE EXTENDED FINITE-ELEMENT METHODComputer Methods in Applied Mechanics and Engineering, 190
N. Sharma, Chad Landis, Pradeep Sharma (2010)
Piezoelectric thin-film superlattices without using piezoelectric materialsJournal of Applied Physics, 108
R. Mindlin, N. Eshel (1968)
On first strain-gradient theories in linear elasticityInternational Journal of Solids and Structures, 4
Sheng Mao, P. Purohit, N. Aravas (2016)
Mixed finite-element formulations in piezoelectricity and flexoelectricityProceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 472
L. Cordes, B. Moran (1996)
Treatment of material discontinuity in the Element-Free Galerkin methodComputer Methods in Applied Mechanics and Engineering, 139
Shuling Hu, S. Shen (2010)
Variational principles and governing equations in nano-dielectrics with the flexoelectric effectScience China Physics, Mechanics and Astronomy, 53
B. Nguyen, X. Zhuang, T. Rabczuk (2019)
NURBS-based formulation for nonlinear electro-gradient elasticity in semiconductorsComputer Methods in Applied Mechanics and Engineering
Feng Deng, Q. Deng, S. Shen (2018)
A Three-Dimensional Mixed Finite Element for FlexoelectricityJournal of Applied Mechanics, 85
M. Majdoub, P. Sharma, T. Çagin (2008)
Enhanced size-dependent piezoelectricity and elasticity in nanostructures due to the flexoelectric effectPhysical Review B, 77
M. Lumentut, I. Howard (2017)
Intrinsic electromechanical dynamic equations for piezoelectric power harvestersActa Mechanica, 228
Wenbin Huang, Xiang Yan, S. Kwon, Shujun Zhang, F. Yuan, Xiaoning Jiang (2013)
Flexoelectric strain gradient detection using Ba 0 . 64 Sr 0 . 36 TiO 3 for sensing
A. Soh, Chen Wanji (2004)
Finite element formulations of strain gradient theory for microstructures and the C0–1 patch testInternational Journal for Numerical Methods in Engineering, 61
J. Yvonnet, Liping Liu (2017)
A numerical framework for modeling flexoelectricity and Maxwell stress in soft dielectrics at finite strainsComputer Methods in Applied Mechanics and Engineering, 313
Xu Liang, Shuling Hu, S. Shen (2017)
Nanoscale mechanical energy harvesting using piezoelectricity and flexoelectricitySmart Materials and Structures, 26
G. Camacho, M. Ortiz (1996)
Computational modelling of impact damage in brittle materialsInternational Journal of Solids and Structures, 33
Xiaojin Fu (2011)
Isogeometric Analysis Toward Integration of CAD and CAEJournal of Shanghai Dianji University
This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution
Bo He, B. Javvaji, X. Zhuang (2019)
Characterizing Flexoelectricity in Composite Material Using the Element-Free Galerkin MethodEnergies
Amir Hosnijeh, Daniel Millán, C. Peco, Marino Balaguer, Irene Vicente (2015)
Revisiting pyramid compression to quantify flexoelectricity: a three-dimensional simulation studyPhysical Review B, 91
P. Fischer, Markus Klassen, J. Mergheim, P. Steinmann, R. Müller (2011)
Isogeometric analysis of 2D gradient elasticityComputational Mechanics, 47
P. Yudin, A. Tagantsev (2013)
Fundamentals of flexoelectricity in solidsNanotechnology, 24
Thanh Nguyen, Sheng Mao, Yao-Wen Yeh, P. Purohit, Michael McAlpine (2013)
Nanoscale FlexoelectricityAdvanced Materials, 25
R. Toupin (1956)
The Elastic DielectricIndiana University Mathematics Journal, 5
Xin-Lin Gao, S. Park (2007)
Variational formulation of a simplified strain gradient elasticity theory and its application to a pressurized thick-walled cylinder problemInternational Journal of Solids and Structures, 44
M. Lumentut, Yi-Chung Shu (2018)
A unified electromechanical finite element dynamic analysis of multiple segmented smart plate energy harvesters: circuit connection patternsActa Mechanica, 229
N. Aravas (2011)
Plane-Strain Problems for a Class of Gradient Elasticity Models—A Stress Function ApproachJournal of Elasticity, 104
Z. Kuang (2008)
Some variational principles in electroelastic media under finite deformationScience in China Series G: Physics, Mechanics and Astronomy, 51
R. Mindlin, D. Gazis (1968)
Polarization Gradient in Elastic DielectricsJournal of Applied Mechanics, 42
S. Nanthakumar, X. Zhuang, Harold Park, T. Rabczuk (2017)
Topology optimization of flexoelectric structuresJournal of The Mechanics and Physics of Solids, 105
Jia-wang Hong, D. Vanderbilt (2013)
First-principles theory and calculation of flexoelectricityPhysical Review B, 88
P. Roy, D. Roy (2019)
Peridynamics model for flexoelectricity and damageApplied Mathematical Modelling
D. Codony, O. Marco, Sonia Fern'andez-M'endez, I. Arias (2019)
An immersed boundary hierarchical B-spline method for flexoelectricityComputer Methods in Applied Mechanics and Engineering
energies Review 1,2 3 3 Xiaoying Zhuang , Binh Huy Nguyen , Subbiah Srivilliputtur Nanthakumar , 3 4 4, Thai Quoc Tran , Naif Alajlan and Timon Rabczuk * Division of Computational Mechanics, Ton Duc Thang University, Ho Chi Minh City 758307, Vietnam; xiaoying.zhuang@tdtu.edu.vn Faculty of Civil Engineering, Ton Duc Thang University, Ho Chi Minh City 758307, Vietnam Institute for Continuum Mechanics, Leibniz Universität Hannover, Appelstr. 11, 30167 Hannover, Germany; binh.nguyen@ikm.uni-hannover.de (B.H.N.); nanthakumar@ikm.uni-hannover.de (S.S.N.); tran@ikm.uni-hannover.de (T.Q.T.) Department of Computer Engineering, College of Computer and Information Sciences, King Saud University, Riyadh 11543, Saudi Arabia; timon.rabczuk@uni-weimar.de * Correspondence: timon.rabczuk@uni-weimar.de Received: 21 January 2020; Accepted: 26 February 2020; Published: 12 March 2020 Abstract: Electromechanical coupling devices have been playing an indispensable role in modern engineering. Particularly, ﬂexoelectricity, an electromechanical coupling effect that involves strain gradients, has shown promising potential for future miniaturized electromechanical coupling devices. Therefore, simulation of ﬂexoelectricity is necessary and inevitable. In this paper, we provide an overview of numerical procedures on modeling ﬂexoelectricity. Speciﬁcally, we summarize a generalized formulation including the electrostatic stress tensor, which can be simpliﬁed to retrieve other formulations from the literature. We further show the weak and discretization forms of the boundary value problem for different numerical methods, including isogeometric analysis and mixed FEM. Several benchmark problems are presented to demonstrate the numerical implementation. The source code for the implementation can be utilized to analyze and develop more complex ﬂexoelectric nano-devices. Keywords: ﬂexoelectricity; numerical methods; modeling 1. Introduction In recent years, huge attention in various research areas, such as material science, mechanical engineering and chemistry, has been paid to ﬂexoelectricity [1,2], an electromechanical coupling phenomenon in all dielectric materials (regardless the material symmetry). It describes the induced electric polarization due to strain gradients. Despite the small coefﬁcients as compared to its counterpart, piezoelectricity, ﬂexoelectric effect has gained more attention and is appreciable thanks to the miniaturized trend, where many electromechanical systems are operated at submicron or nano length scale, leading a tremendous increase in the strain gradient. Due to size effects from strain gradients, ﬂexoelectricity holds promise in micro/nano-electromechanical systems such as sensors [3–5], actuators [6,7] and nano-energy harvesters [8–13]. It can also explain physical phenomena such as strain gradient-driven polarity control [14–16], ﬂexo-photovoltaic effects [17] and ﬂexo-caloric effects [18]. Interested readers are referred to other excellent review articles about ﬂexoelectricity for more details [19–22]. To our best knowledge, the ﬁrst continuum theory for ﬂexoelectricity in dielectric solid was proposed by Maranganti et al. [23], who considered both direct and converse ﬂexoelectric effect in the internal energy that takes strain, strain gradient, polarization and polarization gradient as independent variables. As noted in their work, the formulation can be considered as a complement Energies 2020, 13, 1326; doi:10.3390/en13061326 www.mdpi.com/journal/energies Energies 2020, 13, 1326 2 of 29 for the polarization-gradient theory proposed by Mindlin [24] and later on by Sahin and Host [25]. Based on their formulation, Sharma’s group further studied the enhancement of ﬂexoelectricity in nano-structures [26–30]. It should be remarked that, in the above-mentioned formulation, the electric force, which is the divergence of the Maxwell stress, was not taken into account. The consideration of Maxwell stress may be important in nano-dielectric solid, as pointed out by Hu and Shen [31], who utilized the physical variation principle proposed by Kuang [32–35] so that the electrostatic stress appears naturally from the minimization of the electric enthalpy. The signiﬁcance of the Maxwell stress was also discussed in the uniﬁed energy formulation of continuum magneto-electro-elasticity in the work of Liu [36,37]. Nonetheless, by virtue of the non-locality, i.e., strain gradient or polarization gradient, ﬂexoelectricity in a solid dielectric is governed by a fourth-order partial differential equation. Therefore, modeling ﬂexoelectricity is mostly simpliﬁed to a 1D model or problems with simple geometries with axisymmetry (such as a cylinder or a disk) or plate. As in strain gradient elasticity, the difﬁculty in fully multi-dimensional modeling for ﬂexoelectricity arises from the C continuity in the approximation of the displacement ﬁeld. This issue was solved ﬁrstly by Abdollahi et al. [38] using the local maximum-entropy (LME) meshfree method. Analogous to strain gradient elasticity, mixed FEM was also adopted to solve ﬂexoelectric problems by Mao et al. [39,40]. Isogeometric analysis (IGA) is also a robust approach to ﬂexoelectricity, as in [41–44], thanks to its straightforward handling of higher-order continuity. Alternatively, the Argyris triangular element can also be used as a remedy to the C requirement, as in the work of Yvonnet and Liu [45]. More recently, in a similar spirit to IGA approach, Codony et al. [46] proposed the immersed boundary method on hierarchical B-Spline for ﬂexoelectricity problems with higher efﬁciency. Besides, ﬂexoelectric effect and the associated damage behavior were also approached with peridynamics method by Roy and Roy [47]. As the research ﬁeld of ﬂexoelectricity is rapidly developing, the need of numerical simulations is apparently inevitable. However, there has not been a review on the computational aspects of ﬂexoelectricity yet. Therefore, in this paper, we aim to provide an overview as well as detailed procedure on modeling ﬂexoelectricity. Different numerical methods, including isogeometric analysis, mixed FEM and meshfree, are employed to solve the boundary value problem. The remainder of the paper is structured as follows. In Section 2, motivated by the existence of the Maxwell stress tensor in dielectric solid, we summarize the physical variational principle in which the electrostatic stress tensor appears naturally as the minimization of the electric enthalpy from the work of Hu and Shen [31]. Furthermore, we present an alternative formulation, derived from the electric Gibbs free energy (but with different independent variables), which also results in the electrostatic stress tensor. Noticeably, the two formulations can be considered as a generalized formulation, for which, when the electrostatic stress is omitted, the strong forms from previous literature can be retrieved. Consequently, we obtain the weak forms from Galerkin method for different numerical approaches. Details on the numerical implementation of the discretization forms can be found in Section 3 together with source code wherever possible. Several numerical examples to validate as well as illustrate the numerical implementation are presented in Section 4. 2. Ground Work Formulation 2.1. Physical Variational Principle Hu and Shen [31] extended the variational principle from Toupin [48] by the physical variational principle proposed by Kuang [32–35] to derive the governing equations for ﬂexoelectricity, in which the Maxwell’s stress is taken into account. We brieﬂy summarize the theory and formulation in this section. Assuming the internal energy as a function of strain e , strain gradient e , polarization P ij ij,k i and polarization gradient P , i,j 1 1 1 1 U e , e , P , P = a P P + b P P + c # # + d P # + f # P + h P # + g # # , (1) ij ij,k i i,j ij i j ijkl i,j k,l ijkl ij kl ijk i jk ijkl ij k,l ijkl i jk,l ijklmn ij,k lm,n 2 2 2 2 Energies 2020, 13, 1326 3 of 29 from which the constitutive relations are obtained as ¶U s = = c # + d P + f P (2a) ij ijkl kl ijk k ijkl k,l ¶# ij ¶U t = = h P + g # (2b) ijk ijkl l ijklmn lm,n ¶# ij,k ¶U E = = a P + d # + h # (2c) i ij j ijk jk ijkl jk,l ¶P ¶U V = = b P + f # , (2d) ij ijkl k,l ijkl kl ¶P i,j where s , E , t and V are the (local) stress, electric ﬁeld, higher-order stress and higher-order ij i ijk ij electric ﬁeld, respectively, while a , b , c , d , f , and h are the material tensors [32]. In other ij ijkl ijkl ijk ijkl ijkl words, the variation of internal energy can be written as dU = s d# + t d# + E dP + V dP . (3) ij ij ijk ij,k i i ij i,j The energy density is the sum of the internal energy and the internal of the Maxwell self-ﬁeld [24] W = U + # f f , (4) 0 ,i ,i in which the electric Maxwell self-ﬁeld can be deﬁned as the gradient of the electric potential f MS E = f . (5) ,i Next, we can deﬁne the electric enthalpy as H = W E D = U # f f + f P (6) i i 0 ,i ,i ,i i We consider a dielectric solid occupying a volume v bounded by a surface a that separates from the 0 0 surrounding environment v (considered to be air in this study). Taking v = v + v , the equilibrium of such system in the the isothermal process can be deﬁned by the variational principle d H dv = dW, (7) ext f where W is the external work done by body force f , external electric ﬁeld E and free charge r on the solid body and traction t , double-traction t , surface charge s and higher-order electric ﬁeld v ¯ on i i i the surface boundary. Hence, its variation dW is deﬁned as R R R R R R ext n f dW = f du + E dP dv + t du da + t ¯ D (du ) da + v ¯ dP da s df da r df dv , (8) s t V f i i i i i i i i i v i a a a a v in which the operator D is explained in the following section. Within physical variational principle approach, variation of the electric enthalpy can be further expressed as [32–35]. Z Z Z d H dv = dH dv + H du dv , (9) k,k v v v e 1 where H = E P + V P is the energy deducted from the total energy minus the pure deformation i i ij i,j energy. Moreover, the migratory variation is the essential feature of physical variational principle, Energies 2020, 13, 1326 4 of 29 in the sense that the virtual displacement not only produces the strain variation but also causes electric potential and polarization variations, such as [31] df = d f + d f = d f + f du , (10a) f u f ,i i df = d f + d f = d f + f du , (10b) ,i f ,i u ,i f ,i ,ij j dP = d P + d P = d P + P du , (10c) i P i i P i i,j j dP = d P + d P = d P + P du = d P + P du , (10d) i,j P i,j u i,j P i,j i,jk k P i,j i,kj k where d () and d () are the variations produced by the (local) variation of the electric potential and polarization, respectively, whereas d is the migratory variation caused by the virtual displacement. ( ) Subsequently, by substituting Equation (10) into Equation (9) with the use of Equations (6) and (3) and employing Gauss divergence theorem, through a lengthy but straightforward manipulation, the variation of the electric enthalpy can be re-expressed as [31] Z Z h i ES d H dv = s + t s du dv ij,j ijm,mj i ij,j v v nh i h io ES t t ˜ ˜ + s t + s n + D (n )n n t D n t du da m m ij ijm,m j l j ijm ijm i ij l j + t n n D (du ) da ijm m j i Z Z ext + E V + f + E dP dv + V n dP da i ij,j ,i i ij j i v a Z Z (# f + P ) df dv + (# f + P )n df da , (11) 0 ,i i 0 ,i i i ,i v a t n ˜ ˜ where n is the normal vector. D () = (d n n ) ¶ () and D = n ¶ () are the tangential and i kl k l ,l l ,l normal surface differential operators, respectively, which emerges from the treatment of the variation of displacement gradient du on the boundary (see also [31,49]). Notably, in the above equation, as a i,j ES direct result of the variational process, the generalized electrostatic stress s appears as ij 1 1 ES s = V P + (E P + V P )d D f + # f f + f P d kj k,i k k kl k,l ij j ,i ,k ,k ,k k ij ij 2 2 M M = s + t , (12) ij ijm,m M M with s and t the Maxwell stress and the electrostatic stress corresponding to the strain and ij ijm,m polarization gradients, respectively, deﬁned as 1 1 s = E P d D f + f + # f f + f P d (13a) k k ij j ,j ,i ,k ,k ,k k ij ij 2 2 t = V P + V P . (13b) kj k,j kl k,l ijm,m Upon substituting Equations (11) and (8) into Equation (7) and utilizing the arbitrariness of du , df and dP , we can obtain the governing equations and the boundary conditions with the remark that mechanical displacement and polarization do not exist in the surrounding air environment Energies 2020, 13, 1326 5 of 29 ES s t + s + f = 0, in v (14a) ij ijm,m i ij ,j ES s = 0, in v’ (14b) ij,j L ext V + E f + E = 0, in v (14c) ij,j ,i i i # f + P = r , in v (14d) 0 ,ii i,i f = 0, in v’ (14e) ,ii Dirichlet boundary conditions u = u ¯ , on a (15a) i i f = f, on a (15b) n u,n D (u ) = u n = u ¯ n , on a (15c) i i,l l i,l l P = P , on a (15d) i i Neumann boundary conditions ES t t s ˜ ˜ s t + [[s ]] n + D (n )n n t D n t = t , on a (16a) m m ij ijm,m j l j ijm ijm i ij l j t n n = t , on a (16b) ijm j i n (# [[f ]] + P ) = s , on a (16c) i 0 ,i i V n = v ¯ , on a (16d) ij j i 2.2. Alternative Formulation The above boundary value problem is formulated from the total energy density U whose strain # , strain gradient # , polarization P and polarization gradient P are independent variables. ij ij,k i i,j Alternatively, we can also choose the electric Gibbs free energy (which is identical to the electric enthalpy in isothermal system), whose independent variables are strain # , strain gradient # , electric ij ij,k ﬁeld E and electric ﬁeld gradient E [38,40,50]: i i,j 1 1 1 1 G # , # , E , E = k E E b E E + c # # e E # m (E # # E ) (17) ij ij,k i i,j ij i j ijkl i,j k,l ijkl ij kl ijk i jk ijkl k ij,l ij k,l 2 2 2 2 Consequently, the constitutive relations are given as ¶G 1 s = = c # + m E (18a) ij ijkl kl ijkl k,l ¶# 2 ij ¶G 1 t = = m E (18b) ijk ijkl l ¶# 2 ij,k ¶G 1 D = = k E + m # (18c) i ij j ijkl k,l ¶E 2 ¶G 1 Q = = b E m # , (18d) ij ijkl k,l ijkl kl ¶E 2 i,j where s , D , t and V are the (local) stress, electric displacement, higher-order stress and ij i ijk ij higher-order electric displacement, respectively, while k , b , c , e and m are material ij ijkl ijkl ijk ijkl tensors [38]. These relations also imply that the variation of electric Gibbs free energy can be expressed as dG = s d# + t d# D dE + Q dE , (19) ij ij ijk ij,k i i ij i,j Energies 2020, 13, 1326 6 of 29 which then can be used to determine the equilibriums from variational principle d G dv = dW, (20) where W is again the external work, but instead of having the contribution of the external electric ﬁeld and higher-order electric ﬁeld as in the previous formulation, higher-order surface charges q ¯ are prescribed on the boundary, so that its variation dW reads R R R R R R ext n n f ˜ ˜ dW = f du + E dP dv + t du da + t ¯ D (du ) da + q ¯D (df) da s df da r df dv . (21) s t Q f i i i i i i i v i a a a a v By analogy with Equation (10), migratory variation of electric ﬁeld and electric ﬁeld gradient are obtained as dE = d E + d E = d f + E du = d f + E du , (22a) i f i u i f ,i i,p p f ,i p,i p dE = d E + d E = d E + E du = d E + E du (22b) f u f p f p i,j i,j i,j i,j i,jp i,j i,pj Now, substituting Equation (22) into Equation (20) with the use of Equation (19) and employing Gauss divergence theorem, the physical variational the electric Gibbs free energy can be obtained Z Z M M d G dv = s + t s t du dv ij,j i ijk,kj ij,j ijm,mj v v M M t t ˜ ˜ + s n t n + s n + t n + D (n )n n t D n t du da m m ij j ijk,k j j j l j ijm ijm i ij ijm,m l j Z Z + t n n D (du ) da + D + D df dv ijk k j i i,i ij,ji a v Z Z t t n ˜ ˜ ˜ + D n D n + D (n )n n Q D (Q n ) df da + Q n n D (df) da , (23) i i ij,j i l i j ij i ij j ij i j a a M M where s and t are the Maxwell stress and the higher-order electrostatic stress that are direct ij ijm,m results of the variational process and expressed as s = E D E D d (24a) i j k k ij ij t = E Q E Q d . (24b) ij ijm,m i,k kj k,l kl Note that Equations (24a) and (13a) are identical. Finally, by substituting Equations (23) and (21) into Equation (20), we can obtain the governing equations and boundary condition as ES s t + s + f = 0, in v (25a) ij ijm,m i ij ,j ES s = 0, in v’ (25b) ij,j D Q = r , in v (25c) i ij,j ,i D Q = f = 0, in v’ (25d) i ij,j ,ii ,i Dirichlet boundary conditions u = u ¯ , on a (26a) i i f = f, on a (26b) n u,n D (u ) = u n = u ¯ n , on a (26c) i,l l i,l l n f,n ˜ ¯ D (f) = f n = f n , on a (26d) ,l l ,l l Energies 2020, 13, 1326 7 of 29 Neumann boundary conditions ES t t s ˜ ˜ s t + [[s ]] n + D (n )n n t D n t = t , on a (27a) m m ij ijm,m j l j ijm ijm i ij l j t n n = t ¯ , on a (27b) ijk k j i t t s ˜ ˜ [[D Q ]]n + D (n )n n Q D (Q n ) = s , on a (27c) i ij,j i l i j ij ij j l i Q n n = q ¯, on a (27d) ij i j In the above equations, although it is not trivial, we need to bear in mind that in vacuum displacement and polarization do not exist; hence, the quantity D Q is degenerated to f in i ij,j ,i vacuum, which also causes the jump [[D Q ]] on the boundary a . i ij,j For the sake of clarity, we compare the two formulations in Table 1. Let us denote the strong form derived from the internal energy U and the electric Gibbs energy G in Sections 2.1 and 2.2 as P-formulation and E-formulation, respectively. As compared to the classical (local) theory, in both formulations, strain gradient re results its conjugate variable, double stress t , to be cast in the ijk balance of linear momentum as well as surface traction and surface higher-order traction. In addition to the classical Dirichlet boundary condition, the displacement gradient in the normal direction is also speciﬁed on the boundary. In the same manner, the electric ﬁeld gradient dependence in electric Gibbs energy of E-formulation causes its conjugate variable, higher-order electric displacement Q , to appear ij in the modiﬁcation of Gauss law and surface charge. The electric potential gradient in normal direction is also speciﬁed on the boundary. On the other hand, in P-formulation, when polarization gradient is considered, its conjugate variable V enters the intramolecular force balance equation and higher-order ij surface polarization, while the Gauss law and surface charge condition are unchanged as compared to classical theory. Moreover, as mentioned above, both formulations give identical Maxwell stress tensor s , whereas the higher-order Maxwell stress are given in their corresponding electric independent ij ES variables. It is worth noting that the electrostatic stress tensor s in vacuum is reduced to ij ES M s = s = # E E E E d , (28) 0 i j k k ij ij ij which is the usual Maxwell stress tensor in electromagnetism. Furthermore, we remark that different formulations from the literature can be deduced from ES the generalized ones presented here. For instance, in the absence of electrostatic stress s and the ij polarization gradient P (in P-formulation) or the electric ﬁeld gradient E (in E-formulation), we can i,j i,j retrieve the strong form given by Sharma [30], Mao [39], Abdollahi [38], Deng [40] and Nguyen [43]. We summarize the two reduced formulations in Table 2. It should also be noted that, in P-formulation, when the external electric ﬁeld is also omitted, the intramolecular force balance is reduced to E = f ,i or the usual relation between electric ﬁeld and electric potential E = f can be retrieved and the i ,i strong form of two formulations are identical. To this end, which formulation to be used for numerical computation is simply a matter of choice, as the weak form of both formulations are identical and can be written as R R R R R s du + t du + (# E + P )df dv = f du dv + t du da r df dv s df da (29) 0 s s ij i,j ijm i,jm i i ,i i i i i v v a v a for P-formulation and Z Z Z Z Z s du + t du + D df dv = f du dv + t du da r df dv s df da (30) ij i,j ijm i,jm i ,i i i i i s s v v a v a for E-formulation. Note that Equation (29) will be solved with respect to displacement and electric potential ﬁelds. Hence, the electric polarization needs to be expressed in terms of the electric potential 1 1 1 by utilizing the constitutive relation in Equation (2c), i.e., P = a f a d # a h # . i ,j jkl kl jklm kl,m ij ij ij Energies 2020, 13, 1326 8 of 29 On the other hand, as the electric ﬁeld E = f is an independent variable in E-formulation, it i ,i is straightforward to solve for the displacement and electric potential ﬁelds from Equation (30). Therefore, in the following sections, we employ different numerical methods to solve the weak form in Equation (30) of simpliﬁed E-formulation. When full formulations (Table 1) are chosen, the electric Gibbs free energy is still more favorable. In this case, the expression of the electric polarization in terms of electric potential is done in both domain and surface integrals. In other words, the application of the surface polarization in the P-formulation is equivalent to that of the potential normal derivative in E-formulation. t t t ˜ ˜ ˜ Table 1. Comparison of strong form. L = D (n )n n t D n t , M = D (n )n n Q i l m j ijm m ijm l i j ij l j l D (Q n ). ij j U(e, P,re,rP) G(e, E,re,rE) ES ES Momentum in v s t + s + f = 0 s t + s + f = 0 ij ijm,m i ij i ijk,k ij ij ,j ,j 0 ES ES Maxwell in v s = 0 s = 0 ij,j ij,j Balance f f Gauss law in v # f + P = r D Q = r 0 ,ii i,i i ij,j ,i L ext Intramolecular force [24] V + E f + E = 0 ij,j ,i i i Gauss law in v f = 0 f = 0 ,ii ,ii Displacement on a u = u ¯ u = u ¯ i i i i ¯ ¯ Electric potential on a f = f f = f u,n Dirichlet BCs Displacement normal derivative on a u n = u ¯ n u n = u ¯ n i,l l i,l l i,l l i,l l f,n Potential normal derivative on a f n = f n ,l l ,l l Surface polarization on a P = P i i ES ES Surface traction on a s t + [[s ]] n + L = t s t + [[s ]] n + L = t ij ijm,m j i i ij ijk,k j i i ij ij Surface charge on a n # [[f ]] + P = s [[D D ]]n + M = s i 0 ,i i i ij,j i Neumann BCs Higher-order traction on a t n n = t ¯ t n n = t ¯ ijm m j i ijm m j i Higher-order surface charge on a Q n n = q ¯ ij i j Higher-order electric ﬁeld on a V n = v ij j i t t ˜ ˜ Table 2. Comparison of reduced strong form. L = D (n )n n t D n t , M = m m i l j ijm ijm l j t t ˜ ˜ D (n )n n Q D (Q n ). l i j ij ij j l i U(e, P,re) G(e, E,re) References [30,39] References [38,40,43] ES ES Momentum in v s t + s + f = 0 s t + s + f = 0 ij ijm,m i ij ijk,k i ij ij ,j ,j f f Gauss law in v # f + P = r D Q = r 0 ,ii i,i i ij,j Balance ,i L ext Intramolecular force [24] V + E f + E = 0 ij,j ,i i i Gauss law in v f = 0 f = 0 ,ii ,ii Displacement on a u = u ¯ u = u ¯ i i i i ¯ ¯ Dirichlet BCs Electric potential on a f = f f = f u,n Displacement normal derivative on a u n = u ¯ n u n = u ¯ n i,l l i,l l i,l l i,l l s ES ES ¯ ¯ Surface traction on a s t + [[s ]] n + L = t s t + [[s ]] n + L = t ij ijm,m j i i ij j i i ijk,k ij ij Neumann BCs s Surface charge on a n # [[f ]] + P = s [[D D ]]n + M = s i 0 ,i i i ij,j i Higher-order traction on a t n n = t ¯ t n n = t ¯ m m ijm j i ijm j i 3. Computational Modeling The modeling of ﬂexoelectricity is to solve the boundary value problem described in Table 2 or to minimize the weak forms (Equation (29) or Equation (30)). To do this numerically, the existence of strain gradient entails the chosen numerical methods must ensure at least C -continuity of the second-order derivatives of the displacement ﬁeld by means of higher-order (at least C -continuity) shape functions. In the following, we summarize the numerical methods that have been employed in ﬂexoelectricity, including meshless, mixed FE and IGA. Energies 2020, 13, 1326 9 of 29 3.1. Meshfree Method Some of the initial prominent works on computational evaluation of ﬂexoelectric structures is by adopting a meshfree method using local maximum entropy (LME) approximants, as was done by Abdollahi et al. [38]. The meshfree formulation is validated by comparing the analytical and numerical energy conversion factor in a one-dimensional problem. Inspired by the experimental works of Ma and Cross, the LME based meshfree method is adopted to analyze a two-dimensional pyramid structure in [38] and a three-dimensional pyramid structure in [50]. The converse ﬂexoelectricity that leads to strain due to indentation of atomic force microscopy tip onto a SrTiO sample is experimentally studied and simulated using LME-based meshfree method in [51]. Besides the works on LME, Bo He et al. [52] presented analysis of 2D ﬂexoelectric structures by element free Galerkin method which uses moving least square (MLS) approximants. Nevertheless, more works on ﬂexoelectricity using meshfree methods are still required, as both MLS and LME shape functions lack Kronecker delta property, which leads to additional treatment on imposing Dirichlet boundary conditions. Moreover, in the multi-material problem that possesses material interfaces, some special techniques such as domain partitioning [53], Lagrange multiplier [54] and Jump function [55] are utilized to capture material discontinuities. On the other side, the higher-order Dirichlet boundary condition has not been mentioned in the meshfree approach in ﬂexoelectricity. 3.2. Mixed Finite Element Method The mixed ﬁnite element method offers the advantage of capturing C continuity in the domain by utilizing C ﬁnite elements. The ﬁrst work on utilizing mixed ﬁnite element was by Mao et al. [39] for analyzing ﬂexoelectric plate and crack in periodic ﬂexoelectric structure, in which they proposed a mixed FE element denoted as I9-87. The formulation for ﬂexoelectricity adopted in [39] is based on P-formulation and so the electrical degrees of freedom include polarization in addition to electric potential. The mechanical DOFs comprise displacements and relaxed displacement gradients and kinematic constraints are imposed by Lagrange multipliers. Mixed FE was also implemented for E-formulation by the authors of [40,56,57]. Nanthakumar et al. [56] adopted a nine-noded quadrilateral element with 54 DOFs to analyze and optimize two-dimensional ﬂexoelectric structures. Deng et al. [40] presented two triangular elements, T37 and T45, and two quadrilateral elements, Q47 and Q59. A hexahedral element was presented by Deng et al. [57], who utilized the mixed FE to analyze a micro-pyramid structure. The main advantage of mixed FE is its simplicity in imposing higher-order Dirichlet boundary conditions due to the displacement gradient degrees of freedom; hence, it can be directly implemented in FE commercial software. Moreover, due to its C -approximation, mixed FEM can be used to alleviate the modeling of interface and composite ﬂexoelectric structure. However, the method suffers one major drawback: expensive computational cost due to the many degrees of freedom per element, especially in 3D. 3.3. IGA Approach Due to the ﬂexibility of manipulating continuity of the basis functions, IGA has been used in solving ﬂexoelectricity problems for topology optimization [41], soft dielectric material [42] and semiconductor [44]. In those works, NURBS basis functions are chosen to approximate both geometry and ﬁeld variables (displacement and electric potential ﬁelds). By using knot insertion technique [58,59], the desired continuity order of NURBS basis functions can be obtained, as shown in Figure 1. Moreover, by controlling the continuity, IGA approach can be used to model interface [42], 0 1 in which the geometrical approximation across the interface is C -continuity while C -continuity is required in the domains as usual. Besides, in topology optimization works [41], the voids are assumed to be solid material whose density is much smaller than the host material; however, C continuity is imposed across the void–solid phase. A more general treatment for modeling discontinuity across the Energies 2020, 13, 1326 10 of 29 internal boundary can be adopted from the works in [60–62]. The approach has been applied for porous media, composites, multi-ﬁeld and multi-material problems. In these works, the interface element has been applied to construct C -continuous interpolation along both sides of the redeﬁned interface. The element enables the possibility in enforcing the boundary condition via adding mechanical traction deﬁned by a constitutive relation at the internal discontinuity into the equilibrium equation. 1 1 0.9 0.9 0.8 0.8 0.7 0.7 0.6 0.6 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 ξ ξ (a) (b) Figure 1. (a) C continuity quadratic NURBS basis functions from knot vector X = f0, 0, 0, 1/3, 1/3, 2/3, 2/3, 1, 1, 1g; and (b) C continuity quadratic NURBS basis functions from knot vector X = 0, 0, 0, 1/3, 2/3, 1, 1, 1 f g Perfect interfaces result in continuous displacement and traction ﬁelds along the interface, i.e., [[u ]] = 0 (31a) [[t ]] = 0, (31b) where [[]] indicates the jump. Such conditions of weakly-discontinuity displacement can be implemented with suitable enrichment functions as in XFEM [63,64]. Moreover, in strain gradient elasticity problems, additional continuity conditions associated with higher-order terms such as displacement normal derivative and double traction [65,66] [[u n ]] = 0 (32a) i,l l [[t n n ]] = 0, (32b) ijm m j should be satisﬁed along the interface. Furthermore, the interfaces in microstructure problems are often considered to be imperfect (displacement and/or traction ﬁelds are not continuous along the interface) [67,68]. However, these imperfect interface models are limited to local elasticity and do not consider strain gradients. Although the interface element was applied for mechanical problems, the same concept can be inherited and applied for multi-ﬁeld problems such as piezoelectric and ﬂexoelectric composites. It is worth discussing the imposition of the boundary conditions in IGA approach. While the Dirichlet boundary condition can be imposed by means of least-square method [69] due to the lack of Kronecker-delta property of NURBS basis functions, the displacement gradients conditions can be applied by imposing constraint on displacement degrees of freedom with Lagrange multipliers. This idea is often used in solving strain-gradient elasticity [70] but has not been discussed in ﬂexoelectricity modeling works. Numerical implementation of the displacement gradients conditions is shown in the next section. In the reduced formulations without consider electric ﬁeld or polarization gradient, the electric potential gradient or surface polarization needs not to be imposed. However, for modeling ﬂexoelectric R(ξ) R(ξ) Energies 2020, 13, 1326 11 of 29 transducers, it is necessary to model the electrode on the surface of the structure by imposing the equipotential condition. This can be done by utilizing the Lagrange multiplier, as we show in the next section. Note that the equipotential has not been used in previous works on ﬂexoelectricity using IGA. It should be remarked that recently Codony et al. [46] proposed the immersed boundary method on hierarchical B-spline for ﬂexoelectricity problems with higher efﬁciency in terms of complex domain shapes smf boundary conditions handling. Moreover, the non-local conditions on non-smooth boundary (which we have neglected) was also treated. In Table 3, we present a comparison of different numerical methods on some special features in ﬂexoelectricity modeling, including continuity handling, computational cost (via degrees of freedom per element), the ability to enforce higher-order Dirichlet boundary conditions and material discontinuity approximation. Each numerical method comes with pros and cons and should be considered based on the desired application. For instance, mixed FE is more suitable for a simpliﬁed 2D ﬂexoelectric problems, whereas ﬂexoelectric domains consisting of voids or inclusions can be modeled with meshless methods. 3.4. Implementations In this section, we present the computation procedure of IGA and mixed FE approaches to solve problems governed by ﬂexoelectricity, which can serve as a numerical tool for simulating the nano-energy harvesters, sensors or actuators. For the reasons mentioned above, we consider the simpliﬁed E-formulation in Table 2, which was also the choice of the authors of [38,40], hence reasonably suitable for comparison and validation. Mixed FEM The requirement of C continuous basis functions due to the fourth-order PDEs of ﬂexoelectricity excludes (or at least complicates) the use of classical Lagrange polynomials. too ensure C continuity displacement gradients, y are introduced in the formulation and the kinematic constraint, y = ru, is imposed using Lagrange multipliers. The weak form of mechanical and electrostatic equilibrium u 1 u,n 1 is obtained by ﬁnding the u 2 u = u ¯ on a , u 2 H (v) , y 2 y = y on a , y 2 H (v) , f 2 f 1 2 u 1 f = f on a , f 2 H (v) , l 2 L (v), and similarly fﬁu 2 fﬁu = 0 on a , fﬁu 2 H (v) , dy 2 u,n 1 f 1 2 dy = 0 on a , dy 2 H (v) , df 2 df = 0 on a , df 2 H (v) and dl 2 L (v), such that, Z Z dP = (s : #(du) + t ˆ .h ˆ(dy) D E(df)) dW + l : (dyrdu)dW W W Z Z Z (33) + (yru) : dldW du t dG du b dW = 0 W G W Z Z Z #(du) : C : #(u) dW #(du) : e E(f) dW h ˆ(dy) . h E(f) dW W W W Z Z Z . . . . . . + h ˆ(dy) . g . h ˆ(y)dW E(df) e : #(u) dW E(df) h . h ˆ(y) dW W W W (34) Z Z Z Z E(df) k E(f) dW + dy : l dW rdu : l dW dl : ru dW W W W W Z Z Z + dl : y dW = du t dG + du b dW W G W ¶U where t ˆ = and h ˆ is the third-order tensor related to relaxed displacement gradient y as h ˆ = ¶h ˆ (y + y ). The constraint y = ru is imposed in a weighted residual manner by including jk,i ik,j d (yru) : ldW in Equation (33). W Energies 2020, 13, 1326 12 of 29 Table 3. Comparison of different numerical methods. (a) and (b) indicate quadratic elements in 2D and 3D, respectively. Meshfree IGA Mixed FEM Continuity Depends on kernel functions Knot manipulation Nodal displacement gradient D.o.f Depends on domain of inﬂuence 27 (a) or 108 (b) 54 (Q54, [56]); 37 (T37), 47 (Q47), 45 (T45), 59 (Q59) [40]; 233 (3D) [57] Higher-order Dirichlet BCs - Lagrange-multiplier Direct imposition C Interface modeling Lagrange-multiplier, Enriched functions Multi-patch - Energies 2020, 13, 1326 13 of 29 A nine-noded isoparametric element was proposed by Nanthakumar et al. [56]. The degrees of freedoms are u and u at all nodes, relaxed displacement gradients y , y , y and y at the four 1 2 11 12 21 22 corner nodes and electric potential f at the four corner nodes. In addition to these DOFs, the element has four Lagrange multipliers l , l , l and l at the four corner nodes. The ﬂexoelectric element 11 12 21 22 is shown in Figure 2. The biquadratic shape functions, N , are used to interpolate the displacement DOFs u and u . Bilinear shape functions, N , are used to interpolate the relaxed displacement 1 2 gradients, Lagrange multipliers and electric potential, f. The ﬁnite element approximation is as shown below, h h l u (x) = N u ; y (x) = N y ; i i å å i i=1 i=1,3,5,7 (35) h l h l f (x) = N f ; l (x) = N l å i å i i i i=1,3,5,7 i=1,3,5,7 h h l du (x) = N du ; dy (x) = N dy ; å i å i i i i=1 i=1,3,5,7 (36) h l h l df (x) = N df ; dl (x) = N dl å i å i i i i=1,3,5,7 i=1,3,5,7 7 6 8 4 1 2 Figure 2. A nine-noded ﬂexoelectric element. DOFs at circles are u and u ; at squares are y , y , y , 1 2 11 12 21 y and f; and at triangles are l , l , l , and l . 22 11 12 21 22 Substituting Equations (35) and (36), into Equation (34) yields the following expression in terms of nodal DOFs: 0 1 Z Z Z T T T T T T T T @ A du B CB dW u + du B e B dW f + dy H h B dW f u f f u u W W 0 1 Z Z Z T T T T T T @ A + dy H g HdW y + df B eB dW u + df B h H dW y f f W W (37) Z Z Z T T T T l T l l df B kB dW f du B N dW l + dy N N dW l f yu W W W Z Z Z Z T T T l T l l T q T q dl N B dW u + dl N N dW y = du N t dG + du N b dW yu W W G W N Energies 2020, 13, 1326 14 of 29 2 3 ¶N 0 0 0 ¶x 6 7 ¶N 6 7 0 0 0 2 3 6 ¶x 7 ¶N " # l l 0 6 7 1 ¶N 1 ¶N ¶x ¶N 6 0 0 7 6 ¶N 7 2 ¶x 2 ¶x ¶x 6 7 where B = , B = , H = , B = u 4 5 f l yu ¶y ¶N 6 7 ¶N 0 0 0 q q 6 ¶y 7 ¶N ¶N ¶y 6 7 ¶y ¶x ¶N 6 7 0 0 0 ¶y 4 5 l l 1 ¶N 1 ¶N 0 0 2 ¶y 2 ¶y 2 3 ¶N ¶x 6 l 7 ¶N 6 0 7 ¶x 6 7 ¶N 6 7 4 5 ¶y ¶N ¶y The algebraic equations to be solved to obtain displacement and electric potential in a ﬂexoelectric structure are, 2 32 3 2 3 K 0 K K u F uu uf u ul 6 76 7 6 7 0 K K K y 0 6 76 7 6 7 yy yf yl 6 76 7 = 6 7 (38) 4 54 5 4 5 K K K 0 f 0 fu fy ff K K 0 0 l 0 lu ly K = B CB dW ; uu u T T T K = B e B dW = K ; uf f u fu T T K = B N dW = K ul yu l lu T T K = H hB dW = K yf f fy K = B gB dW yy y l T K = N N dW = K yl l ly K = B kB dW ; ff f Z Z qT q F = N tdG + N b dW G W 3.5. IGA Approach Within the Galerkin method, both the trial and test functions are approximated by NURBS basis functions (e) (e) u = N u ˜, du = N du ˜, (39a) u u (e) (e) ˜ ˜ f = N f, df = N df, (39b) f f in which N and N are the matrices containing the local NURBS basis functions of element e deﬁned as u f 2 3 N 0 0 i h i 6 7 1 2 n N = 0 N 0 , N = N N . . . N , (40a) 4 5 i u u u u u 0 0 N h i N = N N . . . N , (40b) 1 2 n Energies 2020, 13, 1326 15 of 29 where N is the NURBS basis function at control point ith and n is the number of control points of an element. In the above equation, u ˜ , f, du ˜ and du ˜ are the nodal control variables associated with displacement and electric potential of the trial and test functions, respectively, which can be written in matrix form as h i 1 1 1 n n n u ˜ = u ˜ u ˜ u ˜ . . . u ˜ u ˜ u ˜ , (41) x y z x y z h i 1 1 1 n n n du ˜ = du ˜ du ˜ du ˜ . . . du ˜ du ˜ du ˜ , (42) x y z x y z h i 1 2 n ˜ ˜ ˜ ˜ f = f f . . . f , (43) h i 1 2 n ˜ ˜ ˜ ˜ df = df df . . . df (44) T T The approximation of strain ﬁeld # = [# # # 2# 2# 2# ] , electric ﬁeld E = [E E E ] , strain 11 22 33 23 13 12 1 2 3 h i h i T T ¶# ¶# ¶# ¶E ¶E ¶E gradient r# = and electric ﬁeld gradient rE = can be obtained by using ¶x ¶y ¶z ¶x ¶y ¶z differential operators B ,B ,H and H on Equation (39) such that u f u f (e) (e) # = B u ˜, d# = B du ˜, (45a) u u (e) (e) ˜ ˜ E = B f, dE = B df, (45b) f f (e) (e) r# = H u ˜, dr# = H du ˜, (45c) u u (e) (e) ˜ ˜ rE = H f, drE = H df, (45d) f f in which B , B , H and H contain the ﬁrst- and second-order partial derivative of the NURBS basis u f u f functions, respectively, deﬁned as 2 3 ¶N 0 0 ¶x 6 7 ¶N i 2 3 6 0 0 7 1 2 n ¶B ¶B ¶B ¶y u u u 6 7 . . . ¶N 6 7 h i ¶x ¶x ¶x 6 7 0 0 1 2 n 6 7 i ¶B ¶B ¶B ¶z 1 2 n 6 u u u 7 B = , B = , H = (46) 6 7 B B . . . B . . . u u u ¶N ¶N u u u i i 4 5 ¶y ¶y ¶y 6 0 7 ¶z ¶y 1 2 n 6 7 ¶B ¶B ¶B u u u 6 ¶N ¶N 7 . . . i i ¶z ¶z ¶z 4 5 ¶z ¶x ¶N ¶N i i ¶y ¶x 2 3 1 2 n 2 3 ¶B ¶B ¶B f f f ¶N . . . ¶x ¶x ¶x ¶x h i 6 7 1 2 n 6 ¶N 7 6 ¶B ¶B ¶B 7 i 1 2 n i f f f B = , B = B B . . . B , H = (47) 4 5 6 7 f f . . . f f f f ¶y ¶y ¶y ¶y 4 5 ¶N 1 2 n ¶B ¶B ¶B f f f ¶z . . . ¶z ¶z ¶z Note that matrices H and H contain second-order derivative with respect to the physical coordinate u f of the basis functions. The details are presented in Appendix A. Now, based on the discretized ﬁnite elements, the weak form in Equation (30) can be re-expressed as h i R R R R R n T T e T T T T T f ¯ ˜ ˜ de s + dre t dE D dv = du f dv + du ˜ t da + f r dv df s ¯ da . (48) (e) s,(e) f,(e) e=1 v v a v a Finally, by substituting Equations (39) and (45) into Equation (48) and utilizing the arbitrariness of the test functions, a system of linear equations can be obtained as " #" # " # K K u ˜ uu uf = , (49) K K f fu ff f Energies 2020, 13, 1326 16 of 29 where the stiffness matrices and force vectors are computed as n Z K = B cB dv , (50a) uu u (e) e=1 n Z T T K = B eB H mB dv , (50b) uf f f u u (e) e=1 T T T T K = B e B B m H dv , (50c) fu u u f f (e) e=1 n Z K = B kB dv , (50d) ff f (e) e=1 n Z Z T T f = N f dv + N t da , (50e) u u u s,(e) v a e=1 Z Z T f T f = N r dv + N t da . (50f) u f f,(e) v a e=1 Note that the ﬂexoelectric tensor m for cubic symmetry crystal k is given as 2 3 m m m 0 0 0 0 0 0 0 0 m 0 0 0 0 m 0 11 12 12 44 44 6 7 m = (51) 4 0 0 0 0 0 m m m m 0 0 0 0 0 0 m 0 05 44 12 11 12 44 0 0 0 0 m 0 0 0 0 m 0 0 m m m 0 0 0 44 44 12 12 11 The symmetry of the ﬂexoelectric tensor can be found in [71,72]. For cubic crystal, the non-zero components of ﬂexoelectric tensor are m , m and m while other coefﬁcients with cycling 1111 1221 1122 indices are equal (e.g., m = m ) and with odd number of indices are zeros (e.g., m = 0) [73]. 1221 1331 1112 In matrix notation, the coefﬁcients m , m and m are denoted as m , m and m , respectively [38]. iiii jiij iijj 11 12 44 4. Numerical Examples 4.1. Hollow Cylinder In this ﬁrst numerical example, we consider an inﬁnitely long hollow cylinder subject to mechanical displacement and electric potential at the inner and outer surface, as depicted schematically in Figure 3a. Taking advantage of the axisymmetry, a quarter of the hollow cylinder is modeled as shown in Figure 3b, in which symmetry boundary conditions are enforced at x = 0 and x = 0. 1 2 This problem is an oft-used benchmark in gradient elasticity [74,75], and has been extended to ﬂexoelectricity [39,40]. In our work, we adapt both internal energy U from [39] and electric Gibbs energy G from [40] as follows 1 1 1 ˆ ˆ U = l# # + m# # + l l# # + 2m# # + f # P + 2 f # P + aP P (52) jj kk jk jk jj,i kk,i jk,i jk,i 12 jj,i i 44 ij,i j i i 2 2 2 and 1 1 1 G = l# # + m# # + l l# # + 2m# # m # E 2m # E kE E , (53) jj kk jk jk jj,i kk,i jk,i jk,i 12 jj,i i 44 ij,i j i i 2 2 2 where l and m are the two Lamé parameters, l is the length scale, m and m are two non-zero 12 44 ˆ ˆ ﬂexoelectric coefﬁcients, and f and f are two ﬂexoelectric coupling coefﬁcients that related to the 1 2 ˆ ˆ ﬂexoelectric coefﬁcients by the susceptibility c = k/# 1 such that m = c f , m = c f , a and k are 0 12 12 44 44 the reciprocal susceptibility and dielectric permittivity, also related to one another by the susceptibility. These material parameters can be found in Table 4. Furthermore, due to the geometry and boundary Energies 2020, 13, 1326 17 of 29 condition, the displacement as well as the electric potential depend only on the radial r of the hollow cylinder such as 1 r r u (r) = c r + c + c I + c K , (54a) r 1 2 3 1 4 1 r l l 0 0 m du u (r) r r f(r) = c ln r + c + + , (54b) ( ) 5 6 k dr r 2 2 where l = l + , I and K are the modiﬁed Bessel functions of the ﬁrst and second kind, 1 1 (l+m)k respectively. The coefﬁcients c , c , c , c , c , c are determined from the boundary conditions 1 2 3 4 5 6 u = u , u = u , (55a) r r r r r=r r=r i o f = f , f = f , (55b) r r r r r=r r=r i o t = t = 0 (55c) rrr rrr r=r r=r i o For numerical modeling, both formulations in Table 2 are solved with IGA, henceforth denoted as IGA(P) and IGA(E) for the formulation obtained from electric enthalpy and Gibbs energy, respectively. In addition, we also demonstrate results from second formulation with mixed FEM. For the IGA approach, we employ 25 25 second-order NURBS basis functions that have C continuity to model a quarter of the hollow cylinder, as shown in Figure 3b, and impose additional symmetric boundary condition on the left and bottom edges. While the homogeneous and inhomogeneous Dirichlet boundary conditions are trivial to impose [69], the boundary conditions for the displacement gradients require constraining the degrees of freedom associated with the boundary control points and their neighbors. For instance, the symmetry boundary conditions on the bottom edge in Figure 3b read u = 0, (56a) x =0 u = 0, (56b) 1,2 x =0 u = 0, (56c) 2,1 x =0 where Equation (56c) is automatically fulﬁlled as a result of Equation (56a). To impose Equation (56b), the nodal parameters (associated with the displacement in x direction u ) of the two bottom layer adjacent-bottom bottom of control points are identical, i.e., u = u . This can be implemented by utilizing x x Lagrange multipliers. On the other hand, the displacement gradient boundary conditions can be imposed effortlessly for the mixed FEM with the nodal displacement gradient. Table 4. Material parameter for the problem depicted in Figure 3a. Y n l m m k 12 44 9 6 6 6 9 139 10 Pa 0.3 2 10 m 1 10 C/m 1 10 C/m 1 10 F/m Energies 2020, 13, 1326 18 of 29 −5 x 10 Element edge Control Points 1.8 1.6 1.4 1.2 0.8 0.6 0.4 0.2 0 0.5 1 1.5 2 1 −5 x 10 (b) IGA discretization (a) Schematic of a hollow cylinder. Figure 3. Schematic of an inﬁnitely long hollow cylinder, where the inner and outer radii are r = 10 m and r = 20 m, respectively. The inner surface is imposed by radial displacement u = 0.045 m 0 r and grounded, i.e., f = 0 V, whereas the radial displacement u = 0.05 m and electric potential r r i o f = 1 V are applied on the outer surface. Symmetry boundary conditions are enforced, such that u = u = u = 0 and u = u = u = 0. 1 1,2 2,1 1,2 2,1 x =0 x =0 x =0 x =0 x =0 x =0 1 1 1 2 2 2 Figure 4 presents numerical results of displacement and electric potential obtained from three approaches, showing good agreement between them. For a closer look, numerical results of radial displacement u , electric potential f, radial strain # and circumferential strain # along the radial rr R qq direction with angle q = 45 are presented and compared with the analytical solutions in Figure 5. Excellent agreement between the numerical and analytical results can be seen. 4.2. Cantilever Beam In this example, we consider a 2D cantilever beam subjected to a point load on the right-end and having two different electrode conﬁgurations, which serve as a demonstration for a nanogenerator [38]. As depicted in Figure 6, the beam is constrained at the left-end, whereas two different electrode conﬁguration are set up for two electrical boundary conditions. In Type-1 boundary condition, the electrode at the right-end of the beam is grounded, i.e., the electrical potential f is prescribed as zero. In Type-2 condition, an electric potential difference is applied between the top and bottom electrode, e.g., the electric potential on the top is ﬁxed to zero, whereas the bottom is induced to have an unknown equipotential V, which is solved by utilizing Lagrange multiplier. The material parameter in this study is single barium titanate (BaTiO ) crystal adopted from [38] and presented in Table 5. The mechanical deformation can be transform into electrical energy via the electromechanical coupling effect in the nanogenerator. The performance of this energy conversion can be evaluated by the electromechanical coupling factor k given as the ratio between the electrostatic energy E and eff elec the elastic energy E elas E k E dv elec 2 v k = = R , (57) eff e C e dv elas 2 v Consequently, the normalized piezoelectric constant is deﬁned as eff e ¯ = , (58) piezo where k is obtained as k without considering the ﬂexoelectric effect. Physically, the normalized piezo eff piezoelectric constant indicates the enhancement of ﬂexoelectric effect on the energy conversion. 2 Energies 2020, 13, 1326 19 of 29 Figure 7 demonstrates the normalized piezoelectric constant e ¯ as the normalized thickness h = e h/m varies under Type-1 conﬁguration. The numerical results in the simpliﬁed case, where only 31 12 k and m are non-zeros to mimic the 1D problem, are compared to the analytical solution where 33 12 e +12(m /h) k = and excellent agreement is shown, which illustrates the great enhancement of eff k Y ﬂexoelectricity with the decrease in beam thickness. In addition, this enhancement can also be seen in the full 2D model, which is also in good agreement with the work of Abdollahi et al. [38]. u r Phi 5.000e-08 1.000e+00 0.66667 4.8333e-8 0.33333 4.6667e-8 4.500e-08 0.000e+00 (a) IGA(E): Distribution of u (b) IGA(E): Distribution of f -8 Radial displacement x10 [m] Electric potential [V] 4.4674 4.5739 4.6804 4.7869 4.8935 - 0.2971 - 0.0375 0.2220 0.4815 0.7410 1.000 5.000 (d) IGA(P): Distribution of f (c) IGA(P): Distribution of u 0.8 0.6 0.4 0.2 46 0 -0.2 (e) Mixed FEM: Distribution of u (f) Mixed FEM: Distribution of f Figure 4. Distribution of radial displacement and electric potential. Energies 2020, 13, 1326 20 of 29 −8 x 10 5.1 Analytical Analytical IGA−E formulation 0.8 IGA−E formulation IGA− P formulation IGA− P formulation Mixed FEM 4.9 Mixed FEM 0.6 4.8 0.4 4.7 0.2 4.6 4.5 −0.2 4.4 −0.4 1 1.2 1.4 1.6 1.8 2 2.2 1 1.2 1.4 1.6 1.8 2 2.2 r −5 r −5 x 10 x 10 −4 −3 x 10 x 10 Analytical IGA−E formulation 4.5 IGA− P formulation Meshfree 3.5 2 Analytical IGA−E formulation IGA− P formulation Meshfree 2.5 −2 −4 2 1 1.2 1.4 1.6 1.8 2 2.2 1 1.2 1.4 1.6 1.8 2 2.2 r r −5 −5 x 10 x 10 Figure 5. Comparison of the exact and numerical results for radial displacement u , electric potential f, radial strain # and circumferential strain # as functions of radius r. rr ff Element edge Control Points -6 0 0.5 1 1.5 -5 (b) x V (a) Figure 6. (a) Schematic of a hollow cylinder. Two different types of electrical boundary condition are imposed. The upper ﬁgure describes Type-1 boundary where zero electric potential is imposed on the right edge. Type-2 boundary is demonstrated in the lower ﬁgure where zero electric potential is applied on the top edge and an unknown equipotential (due to the electrode) on the bottom edge. (b) IGA discretization (top) and mixed FEM discretization (bottom). Table 5. Material parameters for the problem depicted in Figure 6. Y n e m k k F 31 12 11 33 9 2 6 9 9 6 100 10 Pa 0.37 -4.4 C/m 1 10 C/m 1110 F/m 12.4810 F/m 10010 N ǫ u rr r φφ φ Energies 2020, 13, 1326 21 of 29 Analytical 1D (Flexo+Piezo) IGA−1D (Flexo+Piezo) Analytical 1D (Flexo) IGA−1D (Flexo) IGA− 2D (Flexo) Mixed FEM (Flexo+Piezo) Mixed FEM (Flexo) 0 1 2 3 4 5 6 7 8 Normalized thickness Figure 7. Effective normalized piezoelectric constant obtained with Type-1 boundary condition. We further investigate the difference in the normalized effective piezoelectric constant e with two different electrical boundary conditions and present the numerical results in Figure 8. In both boundary condition setups, as the beam thickness decreases, the normalized effective piezoelectric constant increases, which is expected due to the size effect nature of ﬂexoelectricity. Notably, Type-1 boundary condition seems to be more effective than Type-2. Furthermore, both IGA and mixed FEM approach give the same prediction on the performance of two boundary condition types. For a closer look, we present the electric potential distribution for two electrical boundary conditions computed from IGA and mixed FEM in Figure 9 with the height and length of the beam being h = 1.2468 10 (m) and L = 10h, respectively. The results from the two methods are in agreement. For Type-1 boundary condition, the electric potential difference is highest at the left end of the beam, where strain gradient is largest. As for the Type-2 boundary condition, the induced electric potential at the top electrode is around 9 mV. 2.5 1.5 IGA-Open circuit IGA-Close circuit 1 Mixed FEM-Open circuit Mixed FEM-Close circuit 0.5 0 1 2 3 4 5 6 7 8 Figure 8. Effective piezoelectric constant under different electrical boundary conditions. e¯ Normalized effective piezoelectric constant Energies 2020, 13, 1326 22 of 29 (a) IGA-Type-1 boundary-f (b) IGA-Type-2 boundary-f (c) Mixed FEM-Type-1 boundary-f (d) Mixed FEM-Type-2 boundary-f Figure 9. Electric potential distribution under Type-1 and Type-2 electrical boundary conditions. Additionally, we also carried out a three-dimensional simulation of a cantilever beam under Type-2 electrical boundary condition by using IGA. The 3D beam is assumed to have square cross-section, where the width is equal to the height h. The length of the beam as well as the applied force on the right surface are unchanged as in the two-dimensional case. Figure 10 shows the induced electric potential. A similar response has been observed as in 2D case, where largest electric potential is at the left end. Nevertheless, the electric potential at the top electrode has the value of 3.21 V. phi phi 5.113e+00 5.113e+00 2.7161 2.7161 0.3191 0.3191 -2.078e+00 -2.078e+00 (a) Electric potential (b) Electric potential of the middle surface. Figure 10. Electric potential of a 3D bending ﬂexoelectric beam. We also study the converse ﬂexoelectric effect from the electrical boundary condition. While the mechanical point load is omitted, a uniform electric ﬁeld with magnitude jEj = 8 MV/m is applied by setting the electric potential to be V = 8h (MV) on the bottom electrode. In this case, the beam behaves as an actuator and induces bending curvature, which is shown in Figure 11b. The calculated curvature of the bending beam due to converse ﬂexoelectric effect is in good agreement with both numerical results from [38] and the experimental ones from Bursian and OI [76]. Further investigation on the distribution of electric ﬁeld in the thickness direction of the beam is shown in Figure 11a, which is also in agreement with the result from [38]. Non-uniform distribution of electric ﬁeld can be observed, especially at the region near the top and bottom edge of the beam, which yields high electric ﬁeld gradient and consequently induces deformation. Note that the considered ﬂexoelectric beams do not account for the resistive load [10] or the rectifying circuit [77,78]. In [77], a piezoelectric structure is analyzed along with its circuit connections. In addition, in [78], comparisons between the charge type Hamiltonian and voltage type Hamiltonian are performed to identify the output power and voltage of the piezoelectric energy harvester. Development of a fully combined ﬂexoelectric generator and interface circuit might be of signiﬁcance in the future. Energies 2020, 13, 1326 23 of 29 x 10 8.2 Experiment LME meshfree IGA 7.8 7.6 7.4 7.2 -5 -10 -8 -6 -4 -2 10 10 10 10 10 6.8 0 0.5 1 1.5 2 2.5 −6 y(m) x 10 (b) Curvature due to inverse ﬂexoelectric effect (a) Electric ﬁeld across the thickness Figure 11. Beam under electrical load. 4.3. Truncated Pyramid This section presents the numerical results of a two-dimensional truncated pyramid under compression with different type of boundary conditions, as adopted from [38]. Figure 12 depicts the schematic of a truncated pyramid whose top surface is grounded with zero electric potential, whereas the bottom surface is attached with an electrode that results in the equipotential condition on it. Two mechanical boundary conditions are considered, namely the rigid boundary condition where the bottom surface is constrained in x direction and the top surface is subjected to uniform load F, and the ﬂexible boundary condition where both the bottom and top surfaces are subjected to uniform compression load F. For numerical modeling, the problem is discretized by quadratic C -continuity NURBS basis functions, as shown in Figure 12b as well as by mixed FEM as in Figure 12c. −4 x 10 Element edge Control Points 0 0.5 1 1.5 2 1 −3 x 10 (b) IGA discretization (a) Schematic of a truncated pyramid. (c) Mixed FEM discretization Figure 12. Schematic and IGA discretization of a truncated pyramid. Numerical results of electric potential f and compression strain # for the rigid boundary yy condition from IGA and Mixed FEM are presented in Figure 13 and for ﬂexible boundary condition from IGA are shown in Figure 14. Good agreement can be observed between the two methods as well as with the reference results [38]. In rigid model, our IGA and mixed FEM approaches predict E(V/m) 2 Energies 2020, 13, 1326 24 of 29 the equipotential of bottom electrode to be 2.7 and 2.5 mV, respectively, while the reference result is 2.6 mV [38]. Similarly, in ﬂexible model, the bottom equipotential is 5.7 mV, agreeing with the value of 5.6 mV from [38]. Additionally, a three-dimensional version of the truncated pyramid is also studied. The geometric parameters are kept as in 2D case, so that the top and bottom surfaces now have dimension of a a 1 1 and a a , respectively. The ﬂexible boundary condition is considered, where the applied force 2 2 F is identical from the 2D case. The numerical result from IGA is shown in Figure 15. As can be seen in Figure 15c, the electric potential in 3D model is very similar to that of 2D model. However, the induced potential on the top electrode is now 6.918 V, three orders of magnitude greater than the result predicted from 2D model. strain yy Phi 3.416e-06 2.728e-03 0.0018187 -1.7708e-5 0.00090934 -3.5415e-5 -4.971e-05 -2.330e-10 (a) IGA(E) - f (b) IGA(E) - # yy (c) Mixed FEM - f (d) Mixed FEM - # yy Figure 13. Rigid boundary condition. Phi strain Y 5.730e-03 3.416e-06 0.0038201 -2.3075e-5 0.00191 -4.6151e-5 0 -1.263e-10 -6.581e-05 (a) IGA(E) - f (b) IGA(E) - # yy Figure 14. Flexible boundary condition. Energies 2020, 13, 1326 25 of 29 phi Element edge 6.918e+00 Control Points 4.6121 2.306 1.043e-08 (b) Electric potential f (a) IGA mesh phi 6.918e+00 4.6121 2.306 1.043e-08 (c) Electric potential f Figure 15. 3D pyramid. 5. Discussion Since its discovery, ﬂexoelectricity bears signiﬁcance because of its potential to aid developing novel electromechanical coupling devices. To fully utilize this effect, both experimental and simulation approaches are necessary. In view of simulation, in this paper, we provide a generalized formulation of ﬂexoelectricity as the extension from gradient elasticity. Remarkably, the boundary value problems that involve fourth-order PDE in ﬂexoelectricity, necessitating C continuity, are obtained from electric enthalpy as well as electric Gibbs free energy. To overcome the high-order continuity, we employ isogeometric analysis and mixed ﬁnite element, where detailed implementations are reported. We solved benchmark problems using these methods on 2D and 3D ﬂexoelectric problems. The source codes are provided and can be used to study more complex problems. Author Contributions: Conceptualization, X.Z., B.H.N., S.S.N., T.Q.T., N.A., and T.R.; methodology, X.Z., B.H.N., S.S.N., T.Q.T., N.A., and T.R.; software, B.H.N., S.S.N., and T.Q.T.; validation, B.H.N., S.S.N., and T.Q.T.; formal analysis, B.H.N., S.S.N., and T.Q.T.; investigation, B.H.N., S.S.N., and T.Q.T.; writing, X.Z., B.H.N., S.S.N., T.Q.T., N.A., and T.R.; visualization, B.H.N., S.S.N., and T.Q.T.; and supervision, X.Z. and T.R. All authors have read and agreed to the published version of the manuscript. Funding: This research received no external funding. Acknowledgments: The authors extend their appreciation to the Distinguished Scientist Fellowship Program (DSFP) at King Saud University for funding this work. Conﬂicts of Interest: The authors declare no conﬂict of interest. Energies 2020, 13, 1326 26 of 29 Appendix A. Second-Order Derivative of NURBS Basis Functions The computation of second-order derivative with respect to the physical coordinate of the basis functions can be achieved by ﬁrst noticing that the ﬁrst-order derivative of the basis functions with respect to the parametric coordinate is computed from the Jacobian matrix 2 3 2 3 2 3 ¶N ¶y ¶x ¶z ¶N ¶x ¶x ¶x ¶x ¶x 6 7 6 7 6 7 ¶N ¶y ¶N i ¶x ¶z i 6 7 6 7 = (A1) 4 5 ¶h ¶h ¶h ¶h ¶y 4 5 4 5 ¶N ¶y ¶N i ¶x ¶z i ¶z ¶z ¶z ¶z ¶z Next, by taking the second-order derivative with respect to the parametric coordinate and re-arranging, the second-order derivative with respect to the physical coordinate is determined from 2 3 2 3 2 3 2 3 2 2 2 x, y, z, 2y z 2x z 2x y N N x y z x x x ,x ,x ,x ,x ,x ,x i,xx i,xx ,xx ,xx ,xx 6 7 2 2 2 6 7 6 7 6 7 2 3 6 x, y, z, 2y z 2x z 2x y 7 N N x y z h h ,h ,h ,h ,h ,h ,h i,yy i,hh ,hh ,hh ,hh x 6 7 6 7 6 7 6 7 6 7 6 7 6 7 i,x 2 2 2 6 7 N N x y z x, y, z, 2y z 2x z 2x y 6 i,zz7 6 i,zz7 6 ,zz ,zz ,zz76 7 6 z z z ,z ,z ,z ,z ,z ,z 7 6 7 = 6 7 6 74N 5 . (A2) i,y 6 7 6 7 6 7 6 7 N N x y z x x y y z z y z + y z x z + x z x y + x y i,yz i,hz ,hz ,hz ,hz 6 ,h ,h ,h ,h ,h ,h ,h ,h ,h7 ,z ,z ,z ,z ,z ,z ,z ,z ,z 6 7 6 7 6 7 6 7 i,z 4 5 4 5 4 5 N N x y z 4 x x y y z z y z + y z x z + x z x y + x y 5 i,xz i,xz ,zx ,zx ,zx ,z ,x ,z ,x ,z ,x ,z ,x ,x ,z ,z ,x ,x ,z ,z ,x ,x ,z N N x y z x x y y z z y z + y z x z + x z x y + x y i,xy i,xh ,xh ,xh ,xh ,h ,h ,h ,h ,h ,h ,h ,h ,h ,x ,x ,x ,x ,x ,x ,x ,x ,x References 1. Tolpygo, K. Long wavelength oscillations of diamond-type crystals including long range forces. Sov. Phys.-Solid State 1963, 4, 1297–1305. 2. Kogan, S.M. Piezoelectric effect during inhomogeneous deformation and acoustic scattering of carriers in crystals. Sov. Phys.-Solid State 1964, 5, 2069–2070. 3. Huang, W.; Yan, X.; Kwon, S.R.; Zhang, S.; Yuan, F.G.; Jiang, X. Flexoelectric strain gradient detection using Ba0. 64Sr0. 36TiO3 for sensing. Appl. Phys. Lett. 2012, 101, 252903. [CrossRef] 4. Kwon, S.R.; Huang, W.; Zhang, S.; Yuan, F.G.; Jiang, X. Study on a ﬂexoelectric microphone using barium strontium titanate. J. Micromech. Microeng. 2016, 26, 045001. [CrossRef] 5. Merupo, V.I.; Guiffard, B.; Seveno, R.; Tabellout, M.; Kassiba, A. Flexoelectric response in soft polyurethane ﬁlms and their use for large curvature sensing. J. Appl. Phys. 2017, 122, 144101. [CrossRef] 6. Bhaskar, U.K.; Banerjee, N.; Abdollahi, A.; Wang, Z.; Schlom, D.G.; Rijnders, G.; Catalan, G. A ﬂexoelectric microelectromechanical system on silicon. Nat. Nanotechnol. 2016, 11, 263. [CrossRef] 7. Zhang, S.; Liu, K.; Xu, M.; Shen, S. A curved resonant ﬂexoelectric actuator. Appl. Phys. Lett. 2017, 111, 082904. [CrossRef] 8. Rey, A.D.; Servio, P.; Herrera-Valencia, E. Bioinspired model of mechanical energy harvesting based on ﬂexoelectric membranes. Phys. Rev. E 2013, 87, 022505. [CrossRef] 9. Jiang, X.; Huang, W.; Zhang, S. Flexoelectric nano-generator: Materials, structures and devices. Nano Energy 2013, 2, 1079–1092. [CrossRef] 10. Deng, Q.; Kammoun, M.; Erturk, A.; Sharma, P. Nanoscale ﬂexoelectric energy harvesting. Int. J. Solids Struct. 2014, 51, 3218–3225. [CrossRef] 11. Choi, S.B.; Kim, G.W. Measurement of ﬂexoelectric response in polyvinylidene ﬂuoride ﬁlms for piezoelectric vibration energy harvesters. J. Phys. D Appl. Phys. 2017, 50, 075502. [CrossRef] 12. Liang, X.; Hu, S.; Shen, S. Nanoscale mechanical energy harvesting using piezoelectricity and ﬂexoelectricity. Smart Mater. Struct. 2017, 26, 035050. [CrossRef] 13. Zhu, R.; Wang, Z.; Ma, H.; Yuan, G.; Wang, F.; Cheng, Z.; Kimura, H. Poling-free energy harvesters based on robust self-poled ferroelectric ﬁbers. Nano Energy 2018, 50, 97–105. [CrossRef] 14. Lu, H.; Bark, C.W.; De Los Ojos, D.E.; Alcala, J.; Eom, C.B.; Catalan, G.; Gruverman, A. Mechanical writing of ferroelectric polarization. Science 2012, 336, 59–61. [CrossRef] [PubMed] 15. Lu, H.; Liu, S.; Ye, Z.; Yasui, S.; Funakubo, H.; Rappe, A.M.; Gruverman, A. Asymmetry in mechanical polarization switching. Appl. Phys. Lett. 2017, 110, 222903. [CrossRef] Energies 2020, 13, 1326 27 of 29 16. Guo, R.; Shen, L.; Wang, H.; Lim, Z.; Lu, W.; Yang, P.; Gruverman, A.; Venkatesan, T.; Feng, Y.P.; Chen, J. Tailoring Self-Polarization of BaTiO3 Thin Films by Interface Engineering and Flexoelectric Effect. Adv. Mater. Interfaces 2016, 3, 1600737. [CrossRef] 17. Yang, M.M.; Kim, D.J.; Alexe, M. Flexo-photovoltaic effect. Science 2018, 360, 904–907. [CrossRef] 18. Liu, Y.; Chen, J.; Deng, H.; Hu, G.; Zhu, D.; Dai, N. Anomalous thermoelectricity in strained Bi 2 Te 3 ﬁlms. Sci. Rep. 2016, 6, 32661. [CrossRef] 19. Nguyen, T.D.; Mao, S.; Yeh, Y.W.; Purohit, P.K.; McAlpine, M.C. Nanoscale ﬂexoelectricity. Adv. Mater. 2013, 25, 946–974. [CrossRef] 20. Zubko, P.; Catalan, G.; Tagantsev, A.K. Flexoelectric effect in solids. Annu. Rev. Mater. Res. 2013, 43. [CrossRef] 21. Yudin, P.; Tagantsev, A. Fundamentals of ﬂexoelectricity in solids. Nanotechnology 2013, 24, 432001. [CrossRef] [PubMed] 22. Krichen, S.; Sharma, P. Flexoelectricity: A perspective on an unusual electromechanical coupling. J. Appl. Mech. 2016, 83, 030801. [CrossRef] 23. Maranganti, R.; Sharma, N.; Sharma, P. Electromechanical coupling in nonpiezoelectric materials due to nanoscale nonlocal size effects: Green’s function solutions and embedded inclusions. Phys. Rev. B 2006, 74, 014110. [CrossRef] 24. Mindlin, R.D. Polarization gradient in elastic dielectrics. Int. J. Solids Struct. 1968, 4, 637–642. [CrossRef] 25. Sahin, E.; Dost, S. A strain-gradients theory of elastic dielectrics with spatial dispersion. Int. J. Eng. Sci. 1988, 26, 1231–1245. [CrossRef] 26. Majdoub, M.; Sharma, P.; Çagin, ˘ T. Dramatic enhancement in energy harvesting for a narrow range of dimensions in piezoelectric nanostructures. Phys. Rev. B 2008, 78, 121407. [CrossRef] 27. Majdoub, M.; Sharma, P.; Çagin, ˘ T. Erratum: Dramatic enhancement in energy harvesting for a narrow range of dimensions in piezoelectric nanostructures [Phys. Rev. B 78, 121407 (R)(2008)]. Phys. Rev. B 2009, 79, 159901. [CrossRef] 28. Majdoub, M.; Sharma, P.; Cagin, T. Enhanced size-dependent piezoelectricity and elasticity in nanostructures due to the ﬂexoelectric effect. Phys. Rev. B 2008, 77, 125424. [CrossRef] 29. Majdoub, M.S.; Sharma, P.; Çagin, ˘ T. Erratum: Enhanced size-dependent piezoelectricity and elasticity in nanostructures due to the ﬂexoelectric effect [Phys. Rev. B77, 125424 (2008)]. Phys. Rev. B 2009, 79. [CrossRef] 30. Sharma, N.; Landis, C.; Sharma, P. Piezoelectric thin-ﬁlm superlattices without using piezoelectric materials. J. Appl. Phys. 2010, 108, 024304. [CrossRef] 31. Hu, S.; Shen, S. Variational principles and governing equations in nano-dielectrics with the ﬂexoelectric effect. Sci. China Physics, Mech. Astron. 2010, 53, 1497–1504. [CrossRef] 32. Kuang, Z.B. Some problems in electrostrictive and magnetostrictive materials. Acta Mech. Solida Sin. 2007, 20, 219–227. [CrossRef] 33. Kuang, Z. Some variational principles in electroelastic media under ﬁnite deformation. Sci. China Ser. Physics, Mech. Astron. 2008, 51, 1390–1402. [CrossRef] 34. Kuang, Z.b. Internal energy variational principles and governing equations in electroelastic analysis. Int. J. Solids Struct. 2009, 46, 902–911. [CrossRef] 35. Kuang, Z.B. Theory Electroelasticity; Springer: Berlin/Heidelberg, Germany, 2014. 36. Liu, L. On energy formulations of electrostatics for continuum media. J. Mech. Phys. Solids 2013, 61, 968–990. [CrossRef] 37. Liu, L. An energy formulation of continuum magneto-electro-elasticity with applications. J. Mech. Phys. Solids 2014, 63, 451–480. [CrossRef] 38. Abdollahi, A.; Peco, C.; Millan, D.; Arroyo, M.; Arias, I. Computational evaluation of the ﬂexoelectric effect in dielectric solids. J. Appl. Phys. 2014, 116, 093502. [CrossRef] 39. Mao, S.; Purohit, P.K.; Aravas, N. Mixed ﬁnite-element formulations in piezoelectricity and ﬂexoelectricity. Proc. R. Soc. A 2016, 472, 20150879. [CrossRef] 40. Deng, F.; Deng, Q.; Yu, W.; Shen, S. Mixed Finite Elements for Flexoelectric Solids. J. Appl. Mech. 2017, 84, 081004. [CrossRef] 41. Ghasemi, H.; Park, H.S.; Rabczuk, T. A level-set based IGA formulation for topology optimization of ﬂexoelectric materials. Comput. Methods Appl. Mech. Eng. 2017, 313, 239–258. [CrossRef] Energies 2020, 13, 1326 28 of 29 42. Thai, T.Q.; Rabczuk, T.; Zhuang, X. A large deformation isogeometric approach for ﬂexoelectricity and soft materials. Comput. Methods Appl. Mech. Eng. 2018, 341, 718–739. 43. Nguyen, B.; Zhuang, X.; Rabczuk, T. Numerical model for the characterization of Maxwell-Wagner relaxation in piezoelectric and ﬂexoelectric composite material. Comput. Struct. 2018, 208, 75–91. [CrossRef] 44. Nguyen, B.; Zhuang, X.; Rabczuk, T. NURBS-based formulation for nonlinear electro-gradient elasticity in semiconductors. Comput. Methods Appl. Mech. Eng. 2019, 346, 1074–1095. 45. Yvonnet, J.; Liu, L. A numerical framework for modeling ﬂexoelectricity and Maxwell stress in soft dielectrics at ﬁnite strains. Comput. Methods Appl. Mech. Eng. 2017, 313, 450–482. [CrossRef] 46. Codony, D.; Marco, O.; Fernández-Méndez, S.; Arias, I. An Immersed Boundary Hierarchical B-spline method for ﬂexoelectricity. arXiv 2019, arXiv:1902.02567. 47. Roy, P.; Roy, D. Peridynamics model for ﬂexoelectricity and damage. Appl. Math. Model. 2019, 68, 82–112. [CrossRef] 48. Toupin, R.A. The elastic dielectric. J. Ration. Mech. Anal. 1956, 5, 849–915. [CrossRef] 49. Mindlin, R.D.; Eshel, N. On ﬁrst strain-gradient theories in linear elasticity. Int. J. Solids Struct. 1968, 4, 109–124. [CrossRef] 50. Abdollahi, A.; Millán, D.; Peco, C.; Arroyo, M.; Arias, I. Revisiting pyramid compression to quantify ﬂexoelectricity: A three-dimensional simulation study. Phys. Rev. B 2015, 91, 104103. [CrossRef] 51. Abdollahi, A.; Domingo, N.; Arias, I.; Catalan, G. Converse ﬂexoelectricity yields large piezoresponse force microscopy signals in non-piezoelectric materials. Nat. Commun. 2019, 10, 1266. [CrossRef] 52. He, B.; Javvaji, B.; Zhuang, X. Characterizing Flexoelectricity in Composite Material Using the Element-Free Galerkin Method. Energies 2019, 12, 271. [CrossRef] 53. Liu, G.R. Meshfree Methods: Moving Beyond the Finite Element Method; CRC Press: Boca Raton, FL, USA, 2009. 54. Cordes, L.; Moran, B. Treatment of material discontinuity in the element-free Galerkin method. Comput. Methods Appl. Mech. Eng. 1996, 139, 75–89. [CrossRef] 55. Krongauz, Y.; Belytschko, T. EFG approximation with discontinuous derivatives. Int. J. Numer. Methods Eng. 1998, 41, 1215–1233. [CrossRef] 56. Nanthakumar, S.; Zhuang, X.; Park, H.S.; Rabczuk, T. Topology optimization of ﬂexoelectric structures. J. Mech. Phys. Solids 2017, 105, 217–234. [CrossRef] 57. Deng, F.; Deng, Q.; Shen, S. A three-dimensional mixed ﬁnite element for ﬂexoelectricity. J. Appl. Mech. 2018, 85, 031009. [CrossRef] 58. Les Piegl, W.T. The NURBS Book; Springer: Berlin/Heidelberg, Germany, 1997. 59. Cottrell, J.; Hughes, T.; Bazilevs, Y. Isogeometric Analysis: Toward Integration of CAD and FEA; John Wiley & Sons Inc.: Hoboken, NJ, USA, 2009. 60. Camacho, G.T.; Ortiz, M. Computational modelling of impact damage in brittle materials. Int. J. Solids Struct. 1996, 33, 2899–2938. [CrossRef] 61. Schreﬂer, B.A.; Secchi, S.; Simoni, L. On adaptive reﬁnement techniques in multi-ﬁeld problems including cohesive fracture. Comput. Methods Appl. Mech. Eng. 2006, 195, 444–461. [CrossRef] 62. de Borst, R. Fluid ﬂow in fractured and fracturing porous media: A uniﬁed view. Mech. Res. Commun. 2017, 80, 47–57. [CrossRef] 63. Sukumar, N.; Chopp, D.L.; Moës, N.; Belytschko, T. Modeling holes and inclusions by level sets in the extended ﬁnite-element method. Comput. Methods Appl. Mech. Eng. 2001, 190, 6183–6200. [CrossRef] 64. Sukumar, N.; Huang, Z.; Prévost, J.H.; Suo, Z. Partition of unity enrichment for bimaterial interface cracks. Int. J. Numer. Methods Eng. 2004, 59, 1075–1102. [CrossRef] 65. Shu, J.Y.; King, W.E.; Fleck, N.A. Finite elements for materials with strain gradient effects. Int. J. Numer. Methods Eng. 1999, 44, 373–391. [CrossRef] 66. Soh, A.K.; Wanji, C. Finite element formulations of strain gradient theory for microstructures and the C0–1 patch test. Int. J. Numer. Methods Eng. 2004, 61, 433–454. [CrossRef] 67. Gu, S.; He, Q.C. Interfacial discontinuity relations for coupled multiﬁeld phenomena and their application to the modeling of thin interphases as imperfect interfaces. J. Mech. Phys. Solids 2011, 59, 1413–1426. [CrossRef] 68. Liu, J.T.; He, B.; Gu, S.T.; He, Q.C. Implementation of a physics-based general elastic imperfect interface model in the XFEM and LSM context. Int. J. Numer. Methods Eng. 2018, 115, 1499–1525. [CrossRef] 69. Nguyen, V.P.; Anitescu, C.; Bordas, S.P.; Rabczuk, T. Isogeometric analysis: an overview and computer implementation aspects. Math. Comput. Simul. 2015, 117, 89–116. [CrossRef] Energies 2020, 13, 1326 29 of 29 70. Fischer, P.; Klassen, M.; Mergheim, J.; Steinmann, P.; Müller, R. Isogeometric analysis of 2D gradient elasticity. Comput. Mech. 2011, 47, 325–334. [CrossRef] 71. Shu, L.; Wei, X.; Pang, T.; Yao, X.; Wang, C. Symmetry of ﬂexoelectric coefﬁcients in crystalline medium. J. Appl. Phys. 2011, 110, 104106. [CrossRef] 72. Le Quang, H.; He, Q.C. The number and types of all possible rotational symmetries for ﬂexoelectric tensors. Proc. R. Soc. Lond. A Math. Phys. Eng. Sci. 2011, 467, 2369–2386. 73. Hong, J.; Vanderbilt, D. First-principles theory and calculation of ﬂexoelectricity. Phys. Rev. B 2013, 88, 174107. [CrossRef] 74. Gao, X.L.; Park, S. Variational formulation of a simpliﬁed strain gradient elasticity theory and its application to a pressurized thick-walled cylinder problem. Int. J. Solids Struct. 2007, 44, 7486–7499. [CrossRef] 75. Aravas, N. Plane-strain problems for a class of gradient elasticity models—A stress function approach. J. Elast. 2011, 104, 45–70. [CrossRef] 76. Bursian, E.; OI, Z. Changes in curvature of a ferroelectric ﬁlm due to polarization. Sov. Phys. Solid State, USSR 1968, 10, 1121–1124. 77. Lumentut, M.; Shu, Y. A uniﬁed electromechanical ﬁnite element dynamic analysis of multiple segmented smart plate energy harvesters: circuit connection patterns. Acta Mech. 2018, 229, 4575–4604. [CrossRef] 78. Lumentut, M.; Howard, I. Intrinsic electromechanical dynamic equations for piezoelectric power harvesters. Acta Mech. 2017, 228, 631–650. [CrossRef] Sample Availability: Detailed data and codes are available from the authors by request. c 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
Energies – Multidisciplinary Digital Publishing Institute
Published: Mar 12, 2020
You can share this free article with as many people as you like with the url below! We hope you enjoy this feature!
Read and print from thousands of top scholarly journals.
Already have an account? Log in
Bookmark this article. You can see your Bookmarks on your DeepDyve Library.
To save an article, log in first, or sign up for a DeepDyve account if you don’t already have one.
Copy and paste the desired citation format or use the link below to download a file formatted for EndNote
Access the full text.
Sign up today, get DeepDyve free for 14 days.
All DeepDyve websites use cookies to improve your online experience. They were placed on your computer when you launched this website. You can change your cookie settings through your browser.