An electro-mechanical dynamic model for flexoelectric energy harvesters

Flexoelectricity is a universal electro-mechanical coupling effect that occurs in dielectrics of all symmetric groups and becomes dominant at the micro- and nano-scales. It plays an important role in evaluating micro-electro-mechanical systems (MEMS) such as energy harvesters which convert vibrational energy to electric energy. At finer length scales, micro-inertia effects significantly contribute to the behavior of flexoelectric materials due to the mechanical dispersion. Hence, to properly characterize the vibrational behavior of MEMS, a reliable theoretical approach is required accounting for all possible phenomena that affect the output of the system such as voltage or power density. In this work, we present a consistent (dynamic) model and associated computational framework for flexoelectric structures to study the characteristics of the vibrational behavior of energy harvesters showing the dominance of the flexoelectric effect at micro- and nano-scales. In this context, we quantify the impact of the micro-inertia length scale and the flexoelectric dynamic parameter on both frequency and time responses of energy harvesters.


Introduction
Experimental evidence and observations on thin metal films and wires [1][2][3][4], micro-epoxy beams [5,6], and CFRP composites [7] have indicated significant sizedependent phenomena which become dominant at micro-and nanoscopic levels. In continuum mechanics, this effect is commonly modeled by enriching the conventional continuum theories with some material length scale parameter(s). This in turn has lead to the development of higher-order gradient theories such as the couple stress-based elasticity theory proposed by Mindlin and Tiersten [8], Toupin [9] and Koiter [10]; strain gradient elasticity theory based on the seminal work of Mindlin [11], Eringen and Suhubi [12], Green and Rivlin [13], Kröner [14], Mindlin and Eshel [15] and Germain [16], surface stress elasticity theory developed by Mindlin [17], the surface energy model from Gurtin and Murdoch [18,19] and integral nonlocal theories proposed by Eringen [20]. As these theories pro-vide reasonable explanations and their results agreed well with experimental data accounting for size and surface energy effects, they have been applied successfully for characterizing the mechanical behavior of microand nano-scale structures.
Micro-electro-mechanical systems (MEMS) have wide applications in many areas including the automotive industry, smart wearable devices, sports equipment and energy harvesters. Capturing the electromechanical coupling effect at lower length scales requires an appropriate extension of the strain gradient theories. In classic continuum mechanics, the piezoelectric effect describes the linear relation between the induced electric polarization and the applied mechanical strain. Piezoelectricity exists only in non-centrosymmetric dielectrics. Evidence of 'sizedependent piezoelectricity' was confirmed by Multani and Palkar [21] when they found that decreasing grain sizes of PZT decreases the piezoelectric factor. Other experimental observations of size-dependent piezoelectric have been reported from Mishima et al. [22], Buhlmann et al. [23], Cross [24], Harden et al. [25], Baskaran et al. [26], Catalan et al. [27]. This size dependence is commonly named flexoelectricity [24] and is due to the linear coupling between the electric polarization and strain gradients. Different from piezoelectricity which requires non-centrosymmetric constraint, flexoelectricity is a universal phenomenon that exists in all dielectrics including centrosymmetric materials. Flexoelectricity has a wide range of applications including nanodevices [28,29], sensors and actuators [29][30][31] and energy harvesters [32][33][34][35]. An energy harvester is a self-power autonomous electronic device operated by cultivating available existing power from the surrounding environment such as human movements, solar energy, and wind energy. Similar to conventional piezoelectric energy harvesters, flexoelectric-based energy harvesters convert mechanical energy into electrical energy via the mechanical vibrations of actuators.
While there are numerous contributions about piezoelectric energy harvesters, see Raj et al. [36], Caetano et al. [37], and references therein, there are comparatively few studies on flexoelectricity though they have increased in recent years. Deng et al. [32] proposed a flexoelectric Euler-Bernoulli model to analyze the frequency response functions of energy harvesters. The same method was applied by Faroughi et al. [33] to study the influence of tapering geometric properties on the voltage/power outcomes. A nonlinear flexoelectric beam model for energy harvesters was developed by Wang et al. [35] quantifying the influence of different material parameters such as the gradient index. The material volume ratio coefficient on the generated voltage output of functionally graded flexoelectric energy harvesters was studied by Chu et al. [38]. Majdoub et al. [39] reported that a 5 nm BaTiO 3 beam can sustain almost five times the inhomogeneous strain compared to a corresponding macro beam. Deng et al. [32] showed that decreasing the thickness of a flexoelectric beam made of polyvinylidene fluoride (PVDF) from 3 µm to 0.3 µm leads to a two order increase in the magnitude of the energy conversion efficiency; a similar observation was reported from Wang et al. [35] for nonlinear vibrations of PVDF and from Moura et al. for strontium titanate [40].
Polymer-based piezo-and flexoelectric materials such as PVDF and fiber-reinforced polymers improve the structural bending abilities and therefore promise a higher elastic energy. Flexoelectric beam models have been widely used due to their computational efficiency and simplicity. However, these beam models overestimate the natural frequency. Furthermore, most beam models neglect dynamic contributions. Models, that includes both internal length scale and internal inertia gradient affects, are called dynamically consistent [41][42][43]. The micro-inertia effect are important to describe dispersive wave propagations [43]. For dielectric materials, considering a crystal made of at least two different atoms and under mechanical accelerated motion, the unit cell is distorted, which leads to a polarization wave with an amplitude proportional to the exciting acceleration. This phenomenon is called dynamic flexoelectricity [44][45][46] which is nonzero for dielectric materials constructed by different ions due to the mass differences. Deng et al. [34] studied this influence on the performance of energy harvesters employing an Euler-Bernoulli beam model and a reduced-order method. However, the micro-inertia effect was neglected. Furthermore, beam formulations and reduced-order approaches unfortunately cannot capture the complex underlying physics. For vibration problems such as energy harvesting, there are two main issues: firstly, up to date there is no reliable model considering all possible physical phenomena that can influence the output such as voltage and power density. And secondly, no model has been developed that accurately characterize the vibrational behavior of microstructures.
Thus, we formulate a dynamically consistent model for micro-scale piezo-flexoelectric model by including internal length, internal inertia gradient and flexodynamic effects in the context of nonlinear geometric analysis and develop an IGA (Isogeometric Analysis) formulation to investigate the dynamic/vibrational behaviors of micro-scale energy harvesting structures.
2 Theory of couple stress-based piezo-flexoelectricity

Energy formulation
The internal energy density of an elastic dielectric can be expressed as [47][48][49] where is the Green-Lagrange strain i j = 1 2 (F ik F jk −δ i j ) = 1 2 (u i, j + u j,i + u k,i u k, j ) capturing geometrical nonlinearities, u i is the displacement vector, P i the polarization vector, F i j = ∂ x i ∂ X j the deformation gradient,C i jkl the fourth-order elasticity tensor,Ã i j the dielectric tensor,d i jk the piezoelectric tensor and f i jkl denotes the flexoelectric tensor. For a cubic elastic material tensor, we haveC i jkl =C jikl =C i jlk , C 1111 = C 11 ,C 1122 = C 22 andC 2323 = C 44 denoting three independent nonzero components. Taking advantage of Voigt notation, the elasticity tensor C i j , the dielectric tensor A i j and the piezoelectric material d i j can be, respectively, expressed in matrix form as where α i j is the electric susceptibility. The flexoelectric tensor of a cubic material has three independent parameters and is computed as [50,51].
where δ i j is the Kronecker-Delta and δ i jkl is the fourth-order identity tensor which equals to 1 when i = j = k = l = 1. The theory of couple stress model is based on the idea of replacing the third-order strain gradient i j,k by a traceless second-order tensor κ i j = −e iab ja,b [8,9] and then decomposing it into a symmetric and an antisymmetric part as κ i j = κ s i j +κ a i j (see Appendix A for details). Applying this concept, the energy formulation (27) is expressed as The internal elastic energy is where g 1 and g 2 are two material parameters of the nonlocal elasticity, i.e., g 1 = l 2 1 C 12 and g 2 = l 2 2 C 12 and the flexoelectric energy is

Conservation of energy and weak formulations
According to the energy balance principle for isothermal processes, the rate of change of the internal and kinetic energies are equal to the work done by the body forces, surface tractions, and electric field, i.e., P int + P kin = P ext (7) with where τ i and t i indicate the higher-order traction and traction forces, respectively,b i denotes the body force, v i is the velocity vector. The inertia gradient is included to captures the elastic wave dispersion as shown in [41][42][43]. The dynamic flexoelectric effect describes the response of the polarization to the mechanical acceleration m i j [44][45][46]. The polarization dynamics is also included in the above equation, i.e., the dynamic energy related to the rate of the polarization γ i j [44,46,52]. Finally, it can be shown that the initial boundary value problem reads From the local strong form, the integral formulation of the weak form can be obtain by replacing v i from Eq. (9b) by a test function δu i and integrate over the domain . After integration by parts and using the divergence theorem in Eq. (9b) leads to where the acceleration is denoted asü i = a i . Applying the divergence theorem and integral by part on the micro-inertia gives From Askes and Aifantis [43], we know that the term This acceleration gradient is an additional positive term ensuring a positive kinetic energy density and thus stability [43]. Equation (10) can be rewritten as Gauss's law for an isotropic dielectric material and its corresponding boundary conditions are given as The weak formulation of the electrical problem is obtained by multiplying (13b) with a test function δφ yielding Applying integration by parts we finally obtain and Then, we have In this work, the dynamics due to polarization γ i j is neglected for two reasons: • Different from the micro-inertia length scale related to the higher mechanical dispersion, the dynamics of polarization γ i j is related to the higher electrical dispersion. In some applications such as energy harvesting, the mechanical loading is related to the direct flexoelectric effect rather than converse effect due to an applied electrical field. • There has been no work until now measuring or theoretically evaluating the value of γ .

Time integration and Newton iterative solver
In order to solve the nonlinear dynamic equations (12), (18), we consider equilibrium at time step t = n + 1: We take advantage of the implicit Newmark scheme [53,54] to update the velocity v i and the acceleration where the parameters of the Newmark scheme are chosen as β = 1/4 and γ = 1/2. Similarly, the polarization vector is approximated aṡ Linearization of above equations yields where u and φ are the linearization operators with respect to u n+1 i and φ n+1 i , respectively, and Rayleigh damping, which is a linear combination of the mass and the stiffness matrices, is appliedD and ξ 1 and ξ 2 are two damping ratios. The Newton-Raphson iterative is used to solve the linearized system of equations, i.e., ⎡  Let us consider a cantilever beam subjected to distributed bending force q = 1N/m as illustrated in Fig. 1. The beam is clamped on its left edge with u x = u z = φ = 0 and u x,x = 0 (C 1 continuity condition). The material parameters of polyvinylidene difluoride (PVDF) from Table 1 are adopted; the beam's dimensions are as follows: h=10 −6 m, L=50h. The shear modulus is assumed to be zero for a comparative study between numerical results and analytical solutions which are obtained from beam theory as Figure 2 shows excellent agreement between the analytical solutions and numerical results. In the following numerical example, we will illustrate the limit of the beam model and the need for our proposed twodimensional approach-especially for analyzing the vibration behavior.

Relation of natural frequencies and sample thickness
Many previous works have pointed out the dominance of the flexoelectric effect over the piezoelectric effect at the micro-scale [32,35]. To accurately capture the flexodynamic effect at lower scale, the eigenfrequencies should be evaluated properly. The variation of the first eigenfrequency of the cantilever beam from Fig. 3 with respect to sample thickness is investigated. The material parameters of SrTiO 3 from Table 2 with g 2 = 0 are chosen in order to verify our 2D model with the reduced-order beam formulation from Deng et al. [34] where the ratio L/W/h are fixed to be 20/1/1 for all values of thickness h. The results are illustrated in Fig. 4; for h > 20 nm, the beam and 2D models give nearly same values for f. However, when h ≤ 10 nm, the discrepancy of the two models increase from 23.31 to 33.35 and 44.09% for h = 10 nm, h = 7 nm and h = 5 nm, respectively, suggesting the beam model overestimates the natural frequency.

Linear analysis for frequency response and external electrical circuit condition
If the external mechanical excitation is harmonic, i.e., f u (t) =f u E x p(I ωt) and f φ (t) = 0, the solutions have the form u =ū 0 E x p(I ωt) and φ =φ 0 E x p(I ωt) whereū 0 andφ 0 are obtained from Detailed expressions forK uu ,K uu ,K uu are presented in Appendix C. If the beam operates as a capacitor and is connected to an external electrical resistance, the rate of the average electrical displacement D 3 = − 0 ∂φ ∂ξ 3 + P 3 over the thickness can be expressed as a current flow through the resistor by Gauss law With an input harmonic loading, equation (28) is rewritten as  (29) and the power density is computed as Let us consider the frequency response of the base excitation beam from Fig. 3. The same material parameters are used as in Sect. 4.2. The thickness h = 10 nm, the length L = 20 h, the amplitude of excitation is set equal w 0 = 9.81/ω 2 1 as in [33] with ω 1 = 2π f e 1 . To  where the power is a quadratic function of voltage. Figure 5 shows a direct relation between the flexoelectric parameterf and the electrical outputs. As we can see, higher values off give higher voltages and powers; this indicates the importance of seeking higher flexoelectric effect materials for practical applications [60] in order to improve the performance of energy harvester. The effect of the electrical circuit boundary condition with the variation of electrical resistance R can be found in Fig. 6. Increasing the electrical resistance R from 10 to 500 K makes the outcome voltage increases from 1.37 ×10 −9 µV to 3.24 ×10 −8 µV and this totally agrees with Eq. (29). For each value of R, the peaks of output voltage and power density of frequency response curves are obtained at the natural frequency f e 1 = 1.394 GHz. When connecting with an external electric circuit, the beam plays as a capacitor and the external electrical resistance R contributes to the system admittance jωC + 1/R where C is the beam's capacitor. Hence, it does not exhibit a monotonic behavior of the power density when increasing (or decreasing) the load resistance according to Eqs. (29) and (30). To be precise, increasing the electrical resistance R from 10 to 250 K a b However, when R = 500 K the power is decreased to 0.105 W/mm 3 . In this case, R = 250 K is the optimal value of the electrical resistance. In practical applications, it is important to find the optimal value of R where the power density of energy harvesting structure is maximized. The influence of the flexodynamic effect on the frequency response is illustrated in Fig. 7 by varying the flexodynamic parameter m from 4 × 10 −8 Vs 2 /m 2 to 7×10 −8 Vs 2 /m 2 ; this is studied with three cases of nonlocal elastic parameters g 1 = g 2 = h, g 1 = g 2 = 1.3h and g 1 = g 2 = 1.7h. It can be observed from Fig. 7, the flexodynamic effect has an remarkable impact on the outputs of energy harvester. However, the behavior of the influence depends strongly on the elastic length scales. For the case g 1 = g 2 = h (see Fig. 7a, b), the flexodynamicity diminishes the voltage and power density, i.e., higher values of m cause lower values V and P. For the case g 1 = g 2 = 1.7h, the situation is opposite when the flexodynamic effect enhances the outputs, i.e., higher values of m give higher values V andP (see Fig. 7e, f). When g 1 = g 2 = 1.3h, it is observed that increasing m from 4×10 −8 Vs 2 /m 2 to 6×10 −8 Vs 2 /m 2 leads the voltage decrease from 3.13 to 1.25 × 10 −9 µV. However increasing m from 6 × 10 −8 Vs 2 /m 2 to 7×10 −8 Vs 2 /m 2 makes an increase of the voltage from 1.25 to 1.93 × 10 −9 µV (see Fig. 7c). Similar behavior is observed for the power (see Fig. 7d)  Fig. 8 where l d is varied from h to 30 h. Two values of the nonlocal elastic length scale, i.e., g 1 = g 2 = h and g 1 = g 2 = 2h are included in the study. When g 1 = g 2 = h as illustrated in Fig. 8a, b, the magnitudes of voltage and power density get smaller with larger values of l d . More importantly, we finddifferent from the flexodynamic effect-there is a shifting of the system's natural frequency when varying l d . We note that the inclusion of the inertial-related length scale is necessary in order to capture the dispersive behavior of microstructures. When g 1 = g 2 = 2h (see Fig. 8c, d), it can be found that the first natural frequency f e 1 is increased from 1.39 GHz to 1.72 GHz; however, it also leads to less sensitivity of the outputs with respect to l d . Many studies [61,62] have pointed out that the length scale should be small enough and directly related to the lattice space, it also depends on the correlation properties of the medium. It was shown from Askes et al. [63] that for the wave propagation of pulse loads, taking the ratio h/g 1 = h/g 2 = 1 and h/l d = 1 are appropriate choices.

Nonlinear analysis: shoe-mounted flexoelectric energy harvester
Inspired from the well-known work on curved piezoelectric energy harvesters [64], we study a shoemounted micro-scale energy harvester as illustrated in Fig. 9. The heel strike compressing energy from  Table 1 are adopted. Although PVDF materials belong to the orthorhombic mm2 symmetry, the experimental data from Roh et al. [55] show that at room temperature, the assumption of cubic behavior is still acceptable. A tip mass 4 × 10 −8 kg is attached at the midspan at the lower edge (see Fig. 9), and the gravity forces of this tip mass and the beam are considered as well. The sizes of curved beam are as follows: curvature r = 0.5 µm, beam thickness h = 6r/100, width w = 5r and central angle θ = 2.5π/4. The flexodynamic coefficient is computed theoretically following Tagantsev [44], Yudin and Tagantsev [46] and Kvasov and Tagantsev [52] as m i j = χ i j Open electrical circuit conditions, i.e., R → ∞ are assumed. Figure 10a shows the force pattern imitating the applied force on the foot during the running process with F = 3 Hz and P max = 6.5 µN. The impulse function of stepping process is assumed to have the form of Gaussian distribution, i.e., P(t) = P max e −t 2 /2a 2 where a is the full width at half maximum and we choose a 2 = 0.1. The step frequency in human running is f = 3 Hz [65]. Figure 10b shows the vertical displacement u z at the observed point A; a snap through occurs such that u z jumps from −0.098 to −0.265 µm when P is about 4.71 µN. The full model including both piezoelectric and flexoelectric effects give nearly identical results as the model containing only flexoelectric effects; the piezoelectric model predicts only a slightly smaller value of maximum u z . Distributions of the displacement field of shoe-mounted micro-scale energy harvesters at maximum load of P = 6.5 µN and t = 1.172 s are presented in Fig. 11a, b. The highest difference in the electric potential between the upper and lower layers is observed at the beam's mid-span where the loading force is applied (see Fig. 11c). The output voltage under open electrical circuit condition is shown in Fig. 12a. Before the snap through, in the first cycle, the voltage is positive and reaches a maximum of 0.324 mV at t = 0.113 s. The minimum voltage is −3.389 mV at t = 0.172 s after the snap through. This indicates that a higher voltage level can be obtained when operating the flexoelectric energy harvesting in the regime of finite displacement as the inhomogeneity of applied mechanical strain is enhanced. The minimum value of the voltage from piezoelectric material is −0.129 mV. Obviously, at smaller loads, the piezoelectric effect is nearly zero, larger loads make the strain over thickness inhomogeneous and the voltage related to the piezoelectricity over the thickness is nonzero. However, compared to the flexoelectric effect, the contribution of piezoelectricity is insignificant. Our results indicate that the flexoelectric effect dominates the electro-mechanical coupling effect at the micro-scale. The energy conversion factor is an important quantity for evaluating the efficiency of MEMS devices and computed as where W n+1 F denotes the external work and n+1 E is the electrical energy at the time step t = n + 1. For energy harvesters, h ext = 0 and we obtain The time-dependent energy conversion factor (ECF) in four cycles can be found in Fig. 12b for F = 3 Hz. In the first cycle, when the voltage obtains is maximum at t = 0.113 s, the ECF is also maximum with η = 0.49. At t = 0.13 s, the ECF is minimum, η = 0.054, when snap through is observed, indicating the main contribution of the flexoelectric effect to the electro-mechanical coupling. The voltage response of three different excited frequencies (F = 3 Hz, F = 30 Hz and F = 60 Hz) is presented in Fig. 13. Higher frequencies lead to higher inertia and dynamic effects which yields negative voltage in the first half cycle and increases the voltage magnitude in the remaining half a b c Fig. 13 Voltage response was observed in four cycles for three different exciting frequencies. It is observed that higher frequencies make the inertial effect more profound cycle. The size effect is examined by varying the structural radius r , see Fig. 14. The voltage output is studied for five different values for r . Figure 14a  268 mV for r = 0.5 µm and r = 0.7 µm, respectively. This is because for the cases of r = 0.5 µm and r = 0.7 µm, the structures become more compliant; hence, the snap-through effect happens with a smaller value of applied force. After the snap-through, the deformation is altered from compression to dilation; this leads to a redistribution of the strain and strain gradient. Consequently, the voltage sign and the polarization direction are changed.
For larger values of r, the beam is stiffer; hence, the buckling effect does not occur such that the voltage does not change its sign. The maximum voltage values are 0.975, 0.701 and 0.503 mV for r = 1.0, 1.5, and 2.0 µm, respectively. In the full and the 'pure' flexoelectric models (see Fig. 14a, b), the maximum voltage decreases with increasing r . The situation is different for the piezoelectric model, see Fig. 14c. Larger values of r lead to higher maximum voltage levels confirming that increasing the structural dimensions increases the contribution of piezoelectricity and decreases the contribution of flexoelectricity.

Conclusions
In this work, we have developed a nonlinear dynamically consistent model to characterize the vibrational behavior of electro-mechanical coupling structures at the micro-scale by including size effects, flexoelectricity, and their corresponding dispersive behaviors, i.e., the internal inertia gradient effect and dynamic flexoelectric effect. We also studied the influence of these effects on the output of the flexoelectric energy harvester for both time and frequency responses. It is pointed out that decreasing the structural size not only diminishes the piezoelectric effect and enhances the flexoelectric effect but also increases the dynamic flexoelectric effect and the internal inertia gradient effect. The proposed model is expected to serve as a numerical tool in evaluating and characterizing the realistic dynamic behavior of MEMS and vibrational structures such as actuators and sensors.
For energy harvesting applications, the dynamic flexoelectric effect reduces the output voltage but does not notably influence the natural frequencies. The micro-inertia effect on the other hand shifts the natural frequencies and increases the obtained output voltage. In the finite displacement regime, larger mechanical loads increase the voltage and power density of the energy harvester. From the numerical experiments, it was revealed that significant enhancement in the electric voltage is obtained due to the increase in inhomogeneous mechanical strain suggesting the applications of nonlinear operating MEMS. However, the energy converting factor-when the structures undergo finite displacements-is quite small, so improving the performance of microstructure energy harvester is one of the future studies where some shape optimization or topology optimization methods can be applied to optimize the converting factor.
We also indicated that, at the micro-scale level, flexoelectricity is the dominant coupling effect over its piezoelectric counterpart. When decreasing the structural dimensions, the piezoelectric effect becomes less significant and the flexoelectric effect offers a larger contribution to the total electrical response.
Funding Open Access funding enabled and organized by Projekt DEAL. The authors have not disclosed any funding.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflict of interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.

Appendices Appendix A. Couple stress-based flexoelectricity model
In couple stress theory, the third-order strain gradient i j,k is replaced by a traceless second-order tensor κ i j = −e iab ja,b [8,9]. The strain gradient tensor can be decomposed into a symmetric and an antisymmetric part as κ i j = κ s i j + κ a i j . The antisymmetric tensor κ a i j can be represented by its corresponding dual vector χ i = 1 2 e i jk κ a k j . The strain gradient elastic energy density is given by 1 2 where g 1 and g 2 are two material parameters of the nonlocal elasticity, i.e., g 1 = l 2 1 C 12 and g 2 = l 2 2 C 12 . By using κ i j as the strain gradient metric, the flexoelectric tensor isf i jk = −e jab f ikab . Denotingf =f 1122 −f 1212 2 , we havẽ f i jk P i κ jk =f i jk P i (κ a jk + κ s jk ) =f i jk P i κ a jk = 2f P i χ i (A.3) This means the symmetric part of κ i j does not contribute to the flexoelectric coupling effect. The internal energy density is rewritten as [48,49] ψ( i j , κ s i j , χ i ) = ψ int ( i j , κ s i j , χ i ) NURBS are defined as 3) The displacement u and its test function δu belong to the following spaces: U ∈ {u ∈ H 2 ( )|u = u * on u * andN · ∇u = u * * on u * * },