Detection of void and metallic inclusion in 2D piezoelectric cantilever beam using impedance measurements

Abstract

The aim of current work is to improve the existing inverse methodology of void-detection based on a target impedance curve, leading to quick-prediction of the parameters of single circular void. In this work, mode-shape dependent shifting phenomenon of peaks of impedance curve with change in void location has been analyzed. A number of initial guesses followed by an iterative optimization algorithm based on univariate method has been used to solve the problem. In each iteration starting from each initial guess, the difference between the computationally obtained impedance curve and the target impedance curve has been reduced. This methodology has been extended to detect single circular metallic inclusion in 2D piezoelectric cantilever beam. A good accuracy level was observed for detection of flaw radius and flaw-location along beam-length, but not the precise location along beam-width.

Keywords

piezoelectricity / impedance curve / mode shapes / inverse problem / flaw detection / curve shifting

Introduction

Piezoelectric materials deform under electric field and generate an electric potential when subjected to mechanical load. This phenomenon was discovered in 1880 by French physicists Jacques and Pierre Curie. Piezoelectric materials have Perovskite crystal structure ({\mathrm{A}}^{\mathrm{2}+}{\mathrm{B}}^{\mathrm{4}+}{\mathrm{X}}_{\mathrm{3}}^{\mathrm{2}-}) [¹]. At temperature below Curie point, the off-centering of the central cation, B⁴⁺ is responsible for dipole moment and polarization of the crystal.

Nowadays, piezoelectric materials are widely used as sensors and actuators in various industry because of their ability to transform an electrical signal into a mechanical signal and vice versa. This kind of applications involves high strain and strain rate and fatigue loading which may lead to failure due to crack propagation. Hence, it is necessary to monitor their health in a routine basis. Using sophisticated instruments, small holes and other defects in piezoelectric material can be detected experimentally. In this article, we will present a computational approach based on impedance curve to do the same. We will mostly discuss the effect of change in size and location of single circular void and metallic inclusion in 2D piezoelectric cantilever beam on its impedance curve and subsequently discuss how this observation can help to identify void and inclusion through some example problems. In Ref. [²] the location of resonance and anti-resonance peaks of electrical impedance curve was used to locate the flaw. This idea has been implemented here.

Modeling piezoelectric materials is complex due the coupled nature of the problem involving the interaction between electric and mechanical fields. Disparate scales in the number is just one issue. Furthermore, we are interested in detection cracks and flaws which leads to a moving interface problem which poses major challenges to a computational method. In moving interface problems, the interface commonly needs to be adjusted to the discretization and for this purpose efficient remeshing techniques have been developed, see for instance the contributions from Areias and coworkers [³–¹²] or the contributions in Refs. [¹³,¹⁴]. A very powerful alternative method is the extended finite element method (XFEM) [¹⁵,¹⁶] or modifications of it such as the smoothed XFEM [¹⁷,¹⁸], the phantom node method [¹⁹–²⁴] or certain multiscale methods for fracture [²⁵–³¹]. XFEM allows for modeling moving interfaces without adjusting the mesh to the discretization. Though XFEM has mostly applied to fracture problems, there are also several applications of XFEM related to this topic, i.e., inverse analysis and optimization, also for piezoelectric materials [³²–³⁸]. Similar advantages are observed in peridynamics [³⁹,⁴⁰], dual-horizon peridynamics [⁴¹,⁴²], meshfree methods [⁴³–⁴⁹], cracking particles methods [⁵⁰–⁵⁴] or extended meshfree methods [⁵⁵–⁶⁵], which basically adopt the same concept as in XFEM. Another powerful tool to accurately capture complex geometries accurately is isogeometric analysis (IGA) developed by Hughes and coworkers [⁶⁶,⁶⁷] and applied to numerous problems including thin shell analysis [⁶⁸–⁷¹], image registration [⁷²,⁷³], localization [⁷⁴] and also piezo/flexoelectricity [⁷⁵–⁷⁷] or optimization [⁷⁷,⁷⁸]. Furthermore, IGA can not only be used within the finite element framework but it has been proven a powerful to within adaptive analysis [⁷⁹–⁸¹], collocation methods [⁸²] or boundary elements [⁸³], to name only a few. There have been also contributions of combining the advantages of XFEM and IGA [⁸⁴–⁸⁷]. There are also several contributions employing inverse analysis for determining material parameters, see e.g., the contributions in Refs. [⁸⁸–⁹¹] and references therein. However, none of them exploit the advantages of commercial finite element method (FEM) software.

The manuscript is structured as follows: The next section focuses on the identification of impedance characteristics; this section also included the statement of the boundary value problem and finite element discretization. The flaw detection is discussed subsequently before results and discussions are presented. We finish our manuscript with some concluding remarks.

Identification of impedance characteristics

The electrical impedance of a material depends on the frequency of the input voltage. The dependence of the electrical impedance Z on the voltage frequency is usually presented in the form of impedance curve. In this section, the constitutive relations followed by finite element formulation and its solution with appropriate boundary conditions will be expounded to generate impedance curve of a 2D piezoelectric cantilever beam.

Constitutive equations

The linearized constitutive relations for a piezoelectric material (neglecting the effect of thermal strains) [⁹²] are given by

(1)\mathit{\sigma }={\mathit{c}}^E\mathit{S}-{\mathit{e}}^{\mathrm{T}}\mathit{E},

(2)\mathit{D}=\mathit{e}\mathit{S}+{\mathit{\epsilon }}^S\mathit{E},

where \mathit{\sigma } is stress developed in the material, \mathit{S} is mechanical strain developed in the material, \mathit{E} is applied electric field, \mathit{D} is dielectric displacement, {\mathit{c}}^E is modulus of elasticity under constant electric field, \mathit{e} is piezoelectric coupling tensor, and {\mathit{\epsilon }}^S is permittivity under constant strain.

The detailed derivation of the above equations can be found in Ref. [⁹³], Equations (1) and (2) can be written in matrix form as

(3)\left(\begin{array}{l}\mathit{\sigma }\\ \mathit{D}\end{array}\right)=\left[\begin{array}{l}{\mathit{c}}^E\\ \mathit{e}\end{array}\begin{array}{l}{\mathit{e}}^{\mathrm{T}}\\ -{\mathit{\epsilon }}^S\end{array}\right]\left(\begin{array}{l}\mathit{S}\\ -\mathit{E}\end{array}\right).

In general, piezoelectric materials are transversely isotropic and hence the material property matrix should be considered taking into account of required symmetries. For a three-dimensional transversely isotropic material, the constitutive relation may be written as

(4)\left(\begin{array}{l}{\sigma }_{\mathrm{1}}\\ {\sigma }_{\mathrm{2}}\\ {\sigma }_{\mathrm{3}}\\ {\sigma }_{\mathrm{4}}\\ {\sigma }_{\mathrm{5}}\\ {\sigma }_{\mathrm{6}}\\ D_{\mathrm{1}}\\ D_{\mathrm{2}}\\ D_{\mathrm{3}}\end{array}\right)=\left[\begin{array}{lllllllll}{c}_{\mathrm{11}}^E& c_{\mathrm{12}}^E& c_{\mathrm{13}}^E& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& e_{\mathrm{31}}\\ c_{\mathrm{12}}^E& c_{\mathrm{11}}^E& c_{\mathrm{13}}^E& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& e_{\mathrm{31}}\\ c_{\mathrm{13}}^E& c_{\mathrm{13}}^E& c_{\mathrm{33}}^E& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& e_{\mathrm{33}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& c_{\mathrm{44}}^E& \mathrm{0}& \mathrm{0}& \mathrm{0}& e_{\mathrm{15}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& c_{\mathrm{44}}^E& \mathrm{0}& e_{\mathrm{15}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& c_{\mathrm{66}}^E& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& e_{\mathrm{15}}& \mathrm{0}& -{\epsilon }_{\mathrm{11}}^S& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& e_{\mathrm{15}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -{\epsilon }_{\mathrm{11}}^S& \mathrm{0}\\ e_{\mathrm{31}}& e_{\mathrm{31}}& e_{\mathrm{33}}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -{\epsilon }_{\mathrm{33}}^S\end{array}\right]\left(\begin{array}{l}{S}_{\mathrm{1}}\\ S_{\mathrm{2}}\\ S_{\mathrm{3}}\\ S_{\mathrm{4}}\\ S_{\mathrm{5}}\\ S_{\mathrm{6}}\\ -E_{\mathrm{1}}\\ -E_{\mathrm{2}}\\ -E_{\mathrm{3}}\end{array}\right).

In this paper, we consider only the case of two-dimensional systems and hence Eq. (4) may be further simplified as

(5)\left(\begin{array}{l}{\sigma }_{\mathrm{1}}\\ {\sigma }_{\mathrm{3}}\\ {\sigma }_{\mathrm{5}}\\ D_{\mathrm{1}}\\ D_{\mathrm{3}}\end{array}\right)=\left[\begin{array}{lllll}{c}_{\mathrm{11}}^E& c_{\mathrm{13}}^E& \mathrm{0}& \mathrm{0}& e_{\mathrm{31}}\\ c_{\mathrm{13}}^E& c_{\mathrm{33}}^E& \mathrm{0}& \mathrm{0}& e_{\mathrm{33}}\\ \mathrm{0}& \mathrm{0}& c_{\mathrm{44}}^E& e_{\mathrm{15}}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& e_{\mathrm{15}}& -{\epsilon }_{\mathrm{11}}^S& \mathrm{0}\\ e_{\mathrm{31}}& e_{\mathrm{33}}& \mathrm{0}& \mathrm{0}& -{\epsilon }_{\mathrm{33}}^S\end{array}\right]\left(\begin{array}{l}{S}_{\mathrm{1}}\\ S_{\mathrm{3}}\\ S_{\mathrm{5}}\\ -E_{\mathrm{1}}\\ -E_{\mathrm{3}}\end{array}\right).

Finite element discretization

Finite element method for dynamic analysis in piezoelectric material can be found in Ref. [⁹⁴]. Here, finite element discretization will be derived starting from basic energy equation to keep continuity of the theory in this article.

From Hamilton’s principle for non-conservative system, we can write,

(6){\displaystyle \underset{{{\displaystyle t}}_{\mathrm{1}}}{\overset{{{\displaystyle t}}_{\mathrm{2}}}{\int }}\left[\delta L\left(q_i,{\dot{q}}_i\right)+\delta W\right]}\mathrm{d}t=\mathrm{0},

(7)L\left(q_i,{\dot{q}}_i\right)=T\left({\dot{q}}_i\right)-U(q_i),

where L represents Lagrangian, q_i is generalized coordinate, T is kinetic energy, U is potential energy associated with conservative forces, and \delta W is virtual work.

Now,

(8)\delta W=Q_i\delta q_i,

where Q_i is non-conservative force (or charge), and \delta q_i is virtual displacement (or potential).

(9)\delta L\left(q_i,{\dot{q}}_i\right)={\displaystyle \frac{\partial L}{\partial q_i}}\delta q_i+{\displaystyle \frac{\partial L}{\partial {\dot{q}}_i}}\delta {\dot{q}}_i.

Therefore, from Eq. (6), after substituting Eqs. (7)–(9), we get

(10){\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\displaystyle \frac{\partial T}{\partial {\dot{q}}_i}}\right)+{\displaystyle \frac{\partial U}{\partial q_i}}={\displaystyle \frac{\partial W}{\partial q_i}}.

Now, ‘q’ can be replaced with the independent variables, displacement (u) and potential (\phi) respectively to get two separate equations:

(11){\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\displaystyle \frac{\partial T}{\partial {\dot{u}}_i}}\right)+{\displaystyle \frac{\partial U}{\partial u_i}}={\displaystyle \frac{\partial W}{\partial u_i}},

(12){\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\displaystyle \frac{\partial T}{\partial {\dot{\phi }}_i}}\right)+{\displaystyle \frac{\partial U}{\partial {\phi }_i}}={\displaystyle \frac{\partial W}{\partial {\phi }_i}}.

Now, we introduce shape functions. In our analysis, we consider linear shape functions for both displacement and potential fields.

(13)\{u\}=[N_i^u]\{u_i\}.

Now, the expression for strain is

(14)\{S\}={\displaystyle \frac{\partial \{u\}}{\partial \{x\}}}={\displaystyle \frac{\partial [N_i^u]}{\partial \{x\}}}\{u_i\}\Rightarrow \{S\}=[B_i^u]\{u_i\}.

Similarly,

(15)\phi =[N_i^{\phi }]\{{\phi }_i\},

and electric field,

(16)\{E\}=-[B_i^{\phi }]\{{\phi }_i\}.

The expression for kinetic energy,

(17)T={\displaystyle {\int }_{\mathrm{\Omega }}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\rho {\{\dot{u}\}}^{\mathrm{T}}\{\dot{u}\}\mathrm{d}\mathrm{\Omega }}.

Potential energy,

(18)U={\displaystyle {\int }_{\mathrm{\Omega }}{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}}[{\{S\}}^{\mathrm{T}}\{\sigma \}-{\{E\}}^{\mathrm{T}}\{D\}]\mathrm{d}\mathrm{\Omega }.

Virtual work,

(19)\delta W={\displaystyle {\int }_{{\Omega }_1}{\{\delta u\}}^T}\{f_b\}d\Omega +{\displaystyle {\int }_{{\Gamma }_1}{\{\delta u\}}^T}\{f_s\}d\Gamma +{\{\delta u\}}^T\{f_c\}+{\displaystyle {\int }_{{\Omega }_2}\delta \phi {\bar{\sigma }}_b}d\Omega +{\displaystyle {\int }_{{\Gamma }_2}\delta \phi }{\bar{\sigma }}_{s}d\Gamma +\delta \phi {\bar{\sigma }}_c,

where f_{\mathrm{b}} is the body force over {\mathrm{\Omega }}_{\mathrm{1}}, f_{\mathrm{s}} is the surface force over {\mathrm{\Gamma }}_{\mathrm{1}}, f_{\mathrm{c}} is the concentrated point load, {\overline{\sigma }}_{\mathrm{b}} is the body charge over {\mathrm{\Omega }}_{\mathrm{2}}, {\overline{\sigma }}_{\mathrm{s}} is the surface charge over {\mathrm{\Gamma }}_{\mathrm{2}}, and {\overline{\sigma }}_{\mathrm{c}} is the concentrated point charge.

Equations (11) and (12), after substituting the expressions for kinetic energy, potential energy and virtual work, can be written as

(20)[M]\left\{{\ddot{u}}_i\right\}+[K_{uu}]\{u_i\}+[K_{u\phi }]\{{\phi }_i\}=\{f_i\},

(21)[K_{\phi u}]\{u_i\}+[K_{\phi \phi }]\{{\phi }_i\}=\left\{{\bar{\sigma }}_i\right\},

where,

[M]={\displaystyle {\int }_{\mathrm{\Omega }}\rho {[N_i^u]}^{\mathrm{T}}[N_i^u]\mathrm{d}\mathrm{\Omega }},

[K_{uu}]={\displaystyle {\int }_{\mathrm{\Omega }}{[B_i^u]}^{\mathrm{T}}}[c^E][B_i^u]\mathrm{d}\mathrm{\Omega },

[K_{u\phi }]={\displaystyle {\int }_{\mathrm{\Omega }}{[B_i^u]}^{\mathrm{T}}}[e^{\mathrm{T}}][B_i^{\phi }]\mathrm{d}\mathrm{\Omega },

[K_{\phi u}]={\displaystyle {\int }_{\mathrm{\Omega }}{[B_i^{\phi }]}^{\mathrm{T}}}[e][B_i^u]\mathrm{d}\mathrm{\Omega },

[K_{\phi \phi }]={\displaystyle {\int }_{\mathrm{\Omega }}{[B_i^{\phi }]}^{\mathrm{T}}}[-{\epsilon }^S][B_i^{\phi }]\mathrm{d}\mathrm{\Omega },

\{f_i\}={\displaystyle {\int }_{{\mathrm{\Omega }}_{\mathrm{1}}}{[N_i^u]}^{\mathrm{T}}\{f_b\}}\mathrm{d}\mathrm{\Omega }+{\displaystyle {\int }_{{\Gamma }_{\mathrm{1}}}{[N_i^u]}^{\mathrm{T}}}\{f_s\}d\Gamma +{[N_i^u]}^{\mathrm{T}}\{f_c\},

\left\{{\overline{\sigma }}_i\right\}={\displaystyle {\int }_{{\mathrm{\Omega }}_{\mathrm{2}}}{[N_i^{\phi }]}^{\mathrm{T}}}{\overline{\sigma }}_b\mathrm{d}\mathrm{\Omega }+{\displaystyle {\int }_{{{\displaystyle \Gamma }}_{\mathrm{2}}}{[N_i^{\phi }]}^{\mathrm{T}}{\overline{\sigma }}_s}\mathrm{d}\mathrm{\Gamma }+{[N_i^{\phi }]}^{\mathrm{T}}{\overline{\sigma }}_c,

where \mathrm{\Omega } and \mathrm{\Gamma } denote the domain and the boundary while the indices 1 and 2 indicate the matrix and void/inclusion material, respectively.

Boundary conditions

Recalling the constitutive equations derived in 2D form (see Eq. (5) in Section 2.1) we consider 2D cantilever beam in x-z plane with length along ‘x’-axis. All the deformations are in x-z plane and polarization is in ‘z’-direction as shown in Fig. 1. The ratio of the dimension in the directions, ‘x’ to ‘z’ is 8:1. Boundary conditions are given below:

\mathrm{A}\mathrm{t}x=\mathrm{0},u^k=\mathrm{0},k=\{\mathrm{1},\mathrm{3}\},

\mathrm{A}\mathrm{t}z=\mathrm{0},\phi =\mathrm{0},

\mathrm{A}\mathrm{t}z=h,\phi =e^{i\omega t},

where h is the height of the beam.

Complex material parameters

To account for the material loss due to dielectric dissipation, mechanical relaxation and imperfect energy conversion, it is required to consider complex dielectric constant, elastic constants and piezoelectric coefficients respectively [⁹⁵].

Complex material parameters are generally not provided in manufacturers catalogue. They can be found in various research papers readily available. The material parameters of Pz26 used in our analysis, were taken from [⁹⁶]. They are also given below:

c_{\mathrm{11}}^E=(\mathrm{131}+\mathrm{1.03}\mathrm{i})\mathrm{G}\mathrm{P}\mathrm{a},

c_{\mathrm{13}}^E=(\mathrm{91.8}+\mathrm{0.98}\mathrm{i})\mathrm{G}\mathrm{P}\mathrm{a},

c_{\mathrm{33}}^E=(\mathrm{111}+\mathrm{0.82}\mathrm{i})\mathrm{G}\mathrm{P}\mathrm{a},

c_{\mathrm{44}}^E=(\mathrm{21.2}+\mathrm{0.24}\mathrm{i})\mathrm{G}\mathrm{P}\mathrm{a},

e_{\mathrm{31}}=(-\mathrm{2.38}-\mathrm{0.028}\mathrm{i})\mathrm{C/}{\mathrm{m}}^{\mathrm{2}},

e_{\mathrm{33}}=(\mathrm{14.05}-\mathrm{0.13}\mathrm{i})\mathrm{C/}{\mathrm{m}}^{\mathrm{2}},

e_{\mathrm{15}}=(\mathrm{8.18}-\mathrm{0.102}\mathrm{i})\mathrm{C/}{\mathrm{m}}^{\mathrm{2}},

{\epsilon }_{\mathrm{11}}^S=(\mathrm{812}+\mathrm{8.3}\mathrm{i}){\epsilon }_{\mathrm{0}},

{\epsilon }_{\mathrm{33}}^S=(\mathrm{723}+\mathrm{7.3}\mathrm{l}\mathrm{i}){\epsilon }_{\mathrm{0}},

where, {\epsilon }_{\mathrm{0}} is the vacuum permittivity, 8.8541×10⁻¹² F/m (farads per meter). Density of Pz26 is 7700 kg/m³. The elastic properties of the metallic inclusion have been taken as same as the parent piezoelectric material.

Note: If we directly put all the material parameters in system international unit then the condition number of these stiffness matrices will be very high and the equations cannot be solved directly without pre-conditioning. To avoid this difficulty, we introduce the unit of mass kg′, where, 1 kg′ = 10⁹ kg. The derived units should also be scaled properly.

Impedance curve generation

We re-write Eqs. (20) and (21) with ‘k’ as index instead of ‘i’, from Section 2.2,

(22)[M]\left\{{\ddot{u}}_k\right\}+[K_{uu}]\{u_k\}+[K_{u\phi }]\{{\phi }_k\}=\{f_k\},

(23)[K_{\phi u}]\{{{\displaystyle u}}_k\}+[K_{\phi \phi }]\{{\phi }_k\}=\left\{{\bar{\sigma }}_k\right\}.

In our analysis, we do not apply any load, hence,

(24)\{f_k\}=\mathrm{0},

(25)\left\{{\bar{\sigma }}_k\right\}=\mathrm{0.}

At z = h, we apply unit voltage varying harmonically,\phi =\mathrm{1}{e}^{i\omega t}. We assume harmonic solution in steady state situation,

(26)\{u_k\}=\{U_k\}{e}^{i\omega t},\{{\phi }_k\}=\{{\Phi }_k\}{e}^{i\omega t}.

Substituting these equations in Eqs. (22) and Eq. (23), we get

(27)-{\omega }^{\mathrm{2}}[M]\{U_k\}+[K_{uu}]\{U_k\}+[K_{u\phi }]\{{\Phi }_k\}=\mathrm{0},

(28)[K_{\phi u}]\{U_k\}+[K_{\phi \phi }]\{{\Phi }_k\}=\mathrm{0.}

Equations (27) and (28) can be written in matrix form as

(29)\left[\begin{array}{ll}-{\omega }^{\mathrm{2}}[M]+[K_{uu}]& [K_{u\phi }]\\ [K_{\phi u}]& [K_{\phi \phi }]\end{array}\right]\left(\begin{array}{l}[U_k]\\ \{{\Phi }_k\}\end{array}\right)=\left(\begin{array}{l}\mathrm{0}\\ \mathrm{0}\end{array}\right).

The above set of equations are homogeneous. But we apply unit potential (similar to applying displacement in mechanical problem) on the top surface of the beam, hence we will get ‘reaction charge’ (similar to reaction force) at top and bottom surface which can be written as,

(30){\bar{\sigma }}_k={\widehat{\bar{\sigma }}}_k{e}^{i\omega t}.

Now, current,

(31)I_k={\displaystyle \frac{\mathrm{d}{\overline{\sigma }}_k}{\mathrm{d}t}}=i\omega {\overline{\sigma }}_k.

Taking the amplitude part of both sides,

(32){\widehat{I}}_k=i\omega {\widehat{\bar{\sigma }}}_k.

It should be noted that all the phase differences between load and potential or displacement has been taken care while considering complex material matrix. Hence, from now onwards, only the amplitude of the parameters will be indicated.

We can sum up the currents through all the nodes on top surface to get the total current, across the surface,

(33)I_{\mathrm{total}}={\displaystyle \sum {\widehat{I}}_k}.

Now, impedance,

(34)Z={\displaystyle \frac{\mathrm{potential\; difference\; between\; top\; and\; bottom\; surface}}{\mathrm{current\; through\; top\; surface}}}={\displaystyle \frac{\mathrm{1}}{I_{\mathrm{total}}}}.

Clearly impedance, Z depends on the frequency of the potential applied on the top surface. Hence, if we vary frequency over certain range, we can plot impedance v/s frequency graph which is nothing but impedance curve.

Characteristics of impedance curve

At resonance, the deflection of the beam and the mechanical energy of the system becomes very high. For a beam made of piezoelectric material, energy conversion between mechanical and electrical energy takes place, which develops high amount of surface charge, since the amplitude of the potential has been kept fixed. Now according to Eqs. (32) and (33), the current flowing through the beam will be increased significantly and the impedance of the beam will fall down according to Eq. (34). Since the mechanical and electrical behaviors of a dielectric material are not coupled, we will not get any peak on the impedance curve as shown in Fig. 2

Flawdetection with the help of impedance curve

In the previous section, the forward problem was solved for a given set of parameters and the corresponding impedance curve was obtained. The forward is simple as the system of equations leads to a unique solution. However, the problem of identifying the defect through impedance curve is an inverse problem. The inverse problem, in general, will not have a unique solution and hence one needs to solve a constrained optimization problem to identify the required solution. In this section we will check the uniqueness of the inverse solution for a single void and find a technique to solve the inverse problem.

Mode shapes and their effect on impedance curve

The size and location of circular void can be accurately expressed by three parameters, radius (r), x-coordinate along the length and y-coordinate along the height of the beam (Fig. 3).

In this section, the mode shapes and their significance on shifting of impedance curve will be explained.

From Fig. 4(a), the location of the anti-node can be found at approximately 0.008 m. In Fig. 4(b), it can be observed that the impedance curve shifts to the left direction until the void reaches at x = 0.008 m, and then it starts to shift to the right side. This happens because whenever the void shifts towards the anti-node of a mode-shape, the beam becomes more flexible in that particular mode causing a reduction in the natural frequency corresponding to that mode.

In Fig. 5 we observe a phenomenon similar to what happened in first mode. From Fig. 5(a), the location of two anti-nodes and the node can be found approximately at 0.006, 0.052 and 0.026 m, respectively. In Fig. 5(b), we can see that the impedance curve shifts to the left direction, i.e., the natural frequency decreases until the void reaches the anti-node at x = 0.006 m, and then it starts to shift to the right side till x = 0.026 m and then again it shifts to the left till the anti-node at x = 0.052 m, and then again to the right. The justification of this is the same as the first mode.

Again in Fig. 6, we find similar behavior to what we observed for the first two modes. Since each natural frequency depends upon the relative position of the nodes and antinodes of corresponding mode shape with respect to the flaw, there would surely be some natural frequencies whose values would be different for different ‘x’-coordinate of the hole.

Comparison study on shifting behavior of impedance curve with all three parameters

In this section, we try to understand the behavior of the impedance curve with change in the hole parameters, ‘r’, ‘x’ and ‘y’, each at a time keeping the other two fixed.

It should be noted that, the frequency range was kept fixed for same mode against different void parameters to properly analyze the effect of change in the void parameters on the impedance curve. The observations along with some supportive graphs are listed below.

In Fig. 7, the peaks shift in one direction (either+f or ‒f) monotonically with increment of radius for all the three modes. But the same does not happen with the change in x-coordinate (Fig. 8). This has been explained with the help of mode-shapes in Section 3.1.

From Figs. 7–9, we understand that the impedance curve is much more sensitive with change in parameters, ‘x’ and ‘r’ than with change in ‘y’. The reason is that the change in y-coordinate neither changes the total mass of the beam unlike change in ‘r’, nor greatly affect the stiffness by shifting the void across nodes and anti-nodes unlike the change in x-coordinate. It can be concluded that, with a fixed value of ‘y’, for each pair of ‘r’ and ‘x’, the impedance curve is unique. Hence, instead of finding all the parameters at a time, we first assume a value of ‘y’ and locate small bounds in ‘x’ and ‘r’ surrounding the solution using optimization. Within this very small search-space, we then try to optimize all the three parameters with an accuracy up to certain desired level.

Methodology

Figure 10 describes the flowchart of the solution methodology. First, the forward problem is solved by assuming a size (r) and a location of the void (x, y). The solution of the forward problem results in a calculated impedance curve (G) which is then compared with the target impedance curve (T) to calculate the error. The objective is to minimize this error through an iterative process as depicted in the flow chart.

In general, the target impedance curve would be readily available from experimental data. However, we have not conducted any such experiment. Instead, we assume some random flaw parameters and solve the forward problem to obtain the target impedance curve.

The beam has been modelled as 2D plane problem. 4-noded quadrilateral element with 3 degrees of freedom (displacement and potential) has been used to generate mesh in Abaqus.

Optimization function

In Figs. 7 and 8, we find that, with change of the parameters, not only the value of the impedance at resonance gets changed, but also the curves shift along the frequency-axis and the later phenomenon is more prominent than the first one.

Therefore, we calculated the root mean square (rms)-difference of the target and simulated impedance curve along the frequency axis (d_f) and also along the impedance axis (d_i) and we took a combination of these values with more priority to d_f and used this value as the function for optimization algorithm. Thus, we can write the minimization function as

(35)f=k{d}_{\mathrm{f}}+d_{\mathrm{i}},

where, k is a constant. (kd_f) is chosen to at least 10 times the value of d_i. To obtain the value of k, we need to solve two forward problems with two different sets of parameters.

From Figs. 7–9, it can be interpreted that the influence of the parameters on the shifting of peaks in impedance curves is the minimum in first mode. The shifting of peaks increases with higher modes. Hence, in the optimization problem, we should consider more than one natural frequency for minimizing the error. However, the mode shape dependency of shifting direction of the peaks results in many local minima in the optimization function distribution along the direction of x-coordinate because more number of nodes and anti-nodes get involved (Fig. 11).

From Fig. 4(a), we can see that there is only one anti-node present in the first mode, which divides the beam into two segments indicating a maximum of two local minima in optimization function surface generated only from the first mode. From Fig. 5(a), we observe that there are two anti-nodes and one node present in the second mode dividing the beam into four segments, which implies there are maximum four local minima in the optimization function surface considering only second mode. Similarly, from Fig. 6(a), we can understand that there are maximum six possible local minima if we calculate the optimization function considering only third mode. If there are many local minima present in the optimization function distribution, then the computational cost increases because it demands more number of initial guesses to ensure that all of the minima have been considered. We took the first three modes in our consideration. Hence, there are maximum twelve number of local minima possible.

We get similar optimization function distribution with many local minima along x-coordinate, if we introduce a metallic inclusion inside piezoelectric material instead of a void (Fig. 12).

After observing these phenomena, instead of directly applying multilevel coordinate search (MCS) optimization method [⁹⁷], we employed univariate search algorithm, similar to random walk, alternatively in ‘x’ and ‘r’ direction while keeping ‘y’ constant, to find the small bounded regions of ‘x’ and ‘r’ where the minimum lies.

Results and discussion

In this section, we demonstrate the usability of the methodology described in earlier sections through some example problems. Essentially, we solve an inverse problem of finding the location and size of a hole or inclusion in an elastic-piezoelectric beam by comparing the resultant impedance curve with that of a target impedance curve.

Example problem 1: Circular void

In the following, we set a desired accuracy level of Dx/l = 1.25% and Dr/w = 1% while keeping the y-coordinate fixed at neutral axis of the beam. Then MCS optimization algorithm was applied to try to find all the void parameters, but because of oscillatory nature of the solution of ‘y’ this step was discarded. We took the beam dimension as length= 8 cm, height= 1 cm. We kept the boundary of search as, r = [1 mm, 2 mm], x = [3 mm, 77 mm] and y = [3 mm, 7 mm].

In this problem, the actual location and size of the void or the target value of the void parameters are r_target = 0.0011 m, x_target = 0.062 m, y_target = 0.0065 m. In subsequent paragraphs, the solution steps involving optimization will be illustrated.

To consider all the local minima, we begin the search in r-x domain with 13 initial guesses of different ‘x’ values, but radius is kept fixed at 1.5 mm as shown in Fig. 13.

We evaluate the optimization function-value from Eq. (35) for each of the 13 initial guesses and if we plot those function-values we find 5 local minima in the optimization function curve as shown in Fig. 14.

As stated earlier, one of the local minima in Fig. 14, after applying optimization algorithm, would lead to the global minimum i.e., the solution. We keep the void parameters corresponding to all these 5 local minima and discard the rest as shown in Fig. 15.

Now, we take the best guess i.e., the void having minimum function-value among 5 voids (shown by solid blue circle in Fig. 15) and apply the developed optimization algorithm to improve its function-value in r-x domain. After 36 function-evaluations, we find one possible solution which is the local minimum corresponding to that void (shown by solid blue circle in Fig. 16).

Now, we optimize the function-value of the next best guess (13th guess in Fig. 14) and get another local minimum in r-x domain as shown by solid blue circle in Fig. 17. Since this solution is better than the previous one, we discard the solution found from 7th guess (Figs. 18 and 19). We can proceed further for other three initial guesses and find the global minimum in x-r domain as shown in Fig. 20.

Number of function-evaluations: 155

Target–r = 0.0011 m, x = 0.062 m

Result–r = 0.0011344 m, x = 0.06122 m

Therefore, error in ‘r’ is 0.344% of beam width and error in ‘x’ is 0.975% of beam length which are accurate for desired level of accuracy.

Now, we take a small region surrounding the values of x_result and r_result with full variation of y-coordinate and apply MCS algorithm. We found the result oscillating around y-coordinate. The reason is that we cannot predict the variation of impedance curve when we change the value of y-coordinate since the variation is not significant for first three modes of the impedance curve as shown earlier in Fig. 9. Hence, it is not possible to estimate the precise location of the y-coordinate by this method.

Example problem 2: Circular inclusion

With inclusion, our desired level of accuracy was same as for void (Dx/l = 1.25%, Dr/w = 1%). Here, we searched only in r-x domain while keeping y-coordinate fixed at neutral axis of the beam. We took the beam dimension as half of the previous one, length= 4 cm, height= 5 mm. We kept the boundary of search as, r = [0.5 mm, 1 mm], x = [1.2 mm, 38.8 mm] and y = [1.2 mm, 3.8 mm].

As we have discussed earlier, the shifting nature of the impedance curve of a beam with a circular metallic inclusion is similar to a beam with a circular void. Hence, the maximum number of possible local minima in optimization function distribution would be the same for both the cases, beam with void and beam with inclusion if we consider same number of modes. Therefore, we started with 13 initial guesses (like before) of different ‘x’ values, but radius is kept fixed at 0.75 mm.

In this problem, the actual location and size of the inclusion or the target value of the inclusion parameters are r_target = 0.00098 m, x_target = 0.0062 m, y_target = 0.00323 m.

Figure 21 shows the initial guesses along with the target inclusion.

Figures 22–25 describe how the algorithm is approaching towards solution. Each blue circle is representing one inclusion at a time, hence total number of blue circles (solid and dashed) in a figure indicates total number of solutions found at that step and solid blue denotes the best solution at that step.

Number of function-evaluations: 65

Target – r = 0.00098 m, x = 0.0062 m

Result – r = 0.0009488 m, x = 0.00625 m

Therefore, error in ‘r’ is 0.624% of beam width and error in ‘x’ is 0.125% of beam length which are accurate for desired level of accuracy.

Conclusions

The shifting pattern of impedance curve with change in flaw parameters has been presented with proper justification. Based on this behavior, we have established a computational methodology to detect void and metallic inclusion of simple circular geometry in a piezoelectric material. The forward problem (generation of impedance curve for a given system) and the inverse problem (identifying the void/inclusion by comparing the impedance curve with a target curve) have been solved iteratively to get the parameters of the flaw.

There are several limitations and assumptions made on the methodology established in this article. This technique is unable to find the precise value of ‘y’-coordinate of the hole. We have assumed that there is only one flaw and that is circular and not present at the surface of the beam.

We have avoided random search-based optimization to improve the existing inverse methodology by enhancing the smoothness of optimization function distribution, thus the aim of this article has been fulfilled. In the future we aim to account also for uncertainties as proposed e.g., in Ref. [⁹⁸].

References

[1] Galasso F S. Structure, Properties, and Preparation of Perovskite-Type Compounds. Oxford: Pergamon Press, 1969

[2] Nanthakumar S, Lahmer T, Rabczuk T. Detection of flaws in piezoelectric structures using extended fem. International Journal for Numerical Methods in Engineering, 2013, 96(6): 373–389

[3] Areias P, Rabczuk T. Steiner-point free edge cutting of tetrahedral meshes with applications in fracture. Finite Elements in Analysis and Design, 2017, 132: 27–41

[4] Areias P, Rabczuk T, Msekh M. Phase-field analysis of finite-strain plates and shells including element subdivision. Computer Methods in Applied Mechanics and Engineering, 2016, 312: 322–350

[5] Areias P, Rabczuk T, de Sá J C. A novel two-stage discrete crack method based on the screened Poisson equation and local mesh refinement. Computational Mechanics, 2016, 58(6): 1003–1018

[6] Areias P, Msekh M, Rabczuk T. Damage and fracture algorithm using the screened Poisson equation and local remeshing. Engineering Fracture Mechanics, 2016, 158: 116–143

[7] Areias P, Reinoso J, Camanho P, Rabczuk T. A constitutive-based element-by-element crack propagation algorithm with local mesh refinement. Computational Mechanics, 2015, 56(2): 291–315

[8] Areias P, Rabczuk T, Camanho P. Finite strain fracture of 2D problems with injected anisotropic softening elements. Theoretical and Applied Fracture Mechanics, 2014, 72(1): 50–63

[9] Areias P, Rabczuk T, Dias-da Costa D. Element-wise fracture algorithm based on rotation of edges. Engineering Fracture Mechanics, 2013, 110: 113–137

[10] Areias P, Rabczuk T. Finite strain fracture of plates and shells with configurational forces and edge rotations. International Journal for Numerical Methods in Engineering, 2013, 94(12): 1099–1122

[11] Areias P, Rabczuk T, Camanho P. Initially rigid cohesive laws and fracture based on edge rotations. Computational Mechanics, 2013, 52(4): 931–947

[12] Areias P, Reinoso J, Camanho P, de Sá J C, Rabczuk T. Effective 2D and 3D crack propagation with local mesh refinement and the screened Poisson equation. Engineering Fracture Mechanics, 2018, 189: 339–360

[13] Nguyen-Xuan H, Liu G, Bordas S, Natarajan S, Rabczuk T. An adaptive singular ES-FEM for mechanics problems with singular field of arbitrary order. Computer Methods in Applied Mechanics and Engineering, 2013, 253: 252–273

[14] Zhuang X, Huang R, Liang C, Rabczuk T. A coupled thermo-hydro-mechanical model of jointed hard rock for compressed air energy storage. Mathematical Problems in Engineering, 2014, 179169

[15] Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering, 1999, 46(1): 133–150

[16] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. International Journal for Numerical Methods in Engineering, 1999, 45(5): 601–620

[17] Chen L, Rabczuk T, Bordas S, Liu G, Zeng K, Kerfriden P. Extended finite element method with edge-based strain smoothing (ESM-XFEM) for linear elastic crack growth. Computer Methods in Applied Mechanics and Engineering, 2012, 209–212: 250–265

[18] Bordas S, Rabczuk T, Hung N X, Nguyen V, Natarajan S, Bog T, Quan D, Hiep N. Strain smoothing in fem and XFEM. Computers & Structures, 2010, 88(23–24): 1419–1443

[19] Song J H, Areias P M A, Belytschko T. A method for dynamic crack and shear band propagation with phantom nodes. International Journal for Numerical Methods in Engineering, 2006, 67(6): 868–893

[20] Tabarraei A, Song J H, Waisman H. A two-scale strong discontinuity approach for evolution of shear bands under dynamic impact loads. International Journal for Multiscale Computational Engineering, 2013, 11(6): 543–563

[21] Areias P, Song J, Belytschko T. Analysis of fracture in thin shells by overlapping paired elements. International Journal for Numerical Methods in Engineering, 2006, 195: 5343–5360

[22] Chau-Dinh T, Zi G, Lee P S, Rabczuk T, Song J H. Phantom-node method for shell models with arbitrary cracks. Computers and Structures, 2012, 92–93: 242–246

[23] Rabczuk T, Zi G, Gerstenberger A, Wall W. A new crack tip element for the phantom-node method with arbitrary cohesive cracks. International Journal for Numerical Methods in Engineering, 2008, 75(5): 577–599

[24] Vu-Bac N, Nguyen-Xuan H, Chen L, Lee C, Zi G, Zhuang X, Liu G, Rabczuk T. A phantom-node method with edge-based strain smoothing for linear elastic fracture mechanics. Journal of Applied Mathematics, 2013, 978026

[25] Msekh M, Cuong N, Zi G, Areias P, Zhuang X, Rabczuk T. Fracture properties prediction of clay/epoxy nanocomposites with interphase zones using a phase field model. Engineering Fracture Mechanics, 2018, 188: 287–299

[26] Hamdia K, Silani M, Zhuang X, He P, Rabczuk T. Stochastic analysis of the fracture toughness of polymeric nanoparticle composites using polynomial chaos expansions. International Journal of Fracture, 2017, 206(2): 215–227

[27] Silani M, Talebi H, Hamouda A, Rabczuk T. Nonlocal damage modeling in clay/epoxy nanocomposites using a multiscale approach. Journal of Computational Science, 2016, 15: 18–23

[28] Talebi H, Silani M, Rabczuk T. Concurrent multiscale modeling of three dimensional crack and dislocation propagation. Advances in Engineering Software, 2015, 80(C): 82–92

[29] Talebi H, Silani M, Bordas S, Kerfriden P, Rabczuk T. A computational library for multiscale modeling of material failure. Computational Mechanics, 2014, 53(5): 1047–1071

[30] Budarapu P, Gracie R, Yang S W, Zhuang X, Rabczuk T. Efficient coarse graining in multiscale modeling of fracture. Theoretical and Applied Fracture Mechanics, 2014, 69: 126–143

[31] Budarapu P, Gracie R, Bordas S, Rabczuk T. An adaptive multiscale method for quasi-static crack growth. Computational Mechanics, 2014, 53(6): 1129–1148

[32] Zhang C, Nanthakumar S, Lahmer T, Rabczuk T. Multiple cracks identification for piezoelectric structures. International Journal of Fracture, 2017, 206(2): 151–169

[33] Nanthakumar S, Zhuang X, Park H, Rabczuk T. Topology optimization of flexoelectric structures. Journal of the Mechanics and Physics of Solids, 2017, 105: 217–234

[34] Nanthakumar S, Lahmer T, Zhuang X, Park H, Rabczuk T. Topology optimization of piezoelectric nanostructures. Journal of the Mechanics and Physics of Solids, 2016, 94: 316–335

[35] Nanthakumar S, Valizadeh N, Park H, Rabczuk T. Surface effects on shape and topology optimization of nanostructures. Computational Mechanics, 2015, 56(1): 97–112

[36] Nanthakumar S, Lahmer T, Zhuang X, Zi G, Rabczuk T. Detection of material interfaces using a regularized level set method in piezoelectric structures. Inverse Problems in Science and Engineering, 2016, 24(1): 153–176

[37] Nanthakumar S, Lahmer T, Rabczuk T. Detection of multiple flaws in piezoelectric structures using XFEM and level sets. Computer Methods in Applied Mechanics and Engineering, 2014, 275: 98–112

[38] Nanthakumar S, Lahmer T, Rabczuk T. Detection of flaws in piezoelectric structures using extended fem. International Journal for Numerical Methods in Engineering, 2013, 96(6): 373–389

[39] Gerstle W, Sau N, Silling S. Peridynamic modeling of concrete structures. Nuclear Engineering and Design, 2007, 237(12–13): 1250–1258

[40] Silling S. Reformulation of elasticity theory for discontinuities and long-range forces. Journal of the Mechanics and Physics of Solids, 2000, 48(1): 175–209

[41] Ren H, Zhuang X, Rabczuk T. Dual-horizon peridynamics: A stable solution to varying horizons. Computer Methods in Applied Mechanics and Engineering, 2017, 318: 762–782

[42] Ren H, Zhuang X, Cai Y, Rabczuk T. Dual-horizon peridynamics. International Journal for Numerical Methods in Engineering, 2016, 108(12): 1451–1476

[43] Belytschko T, Lu Y, Gu L. Element-free Galerkin methods. International Journal for Numerical Methods in Engineering, 1994, 37(2): 229–256

[44] Fleming M, Chu Y, Moran B, Belytschko T. Enriched element-free Galerkin methods for crack tip fields. International Journal for Numerical Methods in Engineering, 1997, 40(8): 1483–1504

[45] Belytschko T, Tabbara M. Dynamic fracture using element-free Galerkin methods. International Journal for Numerical Methods in Engineering, 1996, 39(6): 923–938

[46] Belytschko T, Krongauz Y, Organ D, Fleming M, Krysl P. Meshless methods: An overview and recent developments. Computer Methods in Applied Mechanics and Engineering, 1996, 139(1–4): 3–47

[47] Rabczuk T, Belytschko T, Xiao S. Stable particle methods based on Lagrangian kernels. Computer Methods in Applied Mechanics and Engineering, 2004, 193(12–14): 1035–1063

[48] Nguyen V, Rabczuk T, Bordas S, Duflot M. Meshless methods: A review and computer implementation aspects. Mathematics and Computers in Simulation, 2008, 79(3): 763–813

[49] Zhuang X, Augarde C, Mathisen K. Fracture modeling using meshless methods and level sets in 3D: Framework and modeling. International Journal for Numerical Methods in Engineering, 2012, 92(11): 969–998

[50] Rabczuk T, Belytschko T. Cracking particles: A simplified meshfree method for arbitrary evolving cracks. International Journal for Numerical Methods in Engineering, 2004, 61(13): 2316–2343

[51] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A simple and robust three-dimensional cracking-particle method without enrichment. Computer Methods in Applied Mechanics and Engineering, 2010, 199(37–40): 2437–2455

[52] Rabczuk T, Samaniego E. Discontinuous modelling of shear bands using adaptive meshfree methods. Computer Methods in Applied Mechanics and Engineering, 2008, 197(6–8): 641–658

[53] Rabczuk T, Areias P, Belytschko T. A simplified mesh-free method for shear bands with cohesive surfaces. International Journal for Numerical Methods in Engineering, 2007, 69(5): 993–1021

[54] Rabczuk T, Belytschko T. A three-dimensional large deformation meshfree method for arbitrary evolving cracks. Computer Methods in Applied Mechanics and Engineering, 2007, 196(29–30): 2777–2799

[55] Bordas S, Rabczuk T, Zi G. Three-dimensional crack initiation, propagation, branching and junction in non-linear materials by an extended meshfree method without asymptotic enrichment. Engineering Fracture Mechanics, 2008, 75(5): 943–960

[56] Rabczuk T, Areias P, Belytschko T. A meshfree thin shell method for nonlinear dynamic fracture. International Journal for Numerical Methods in Engineering, 2007, 72(5): 524–548

[57] Rabczuk T, Bordas S, Zi G. A three-dimensional meshfree method for continuous multiple-crack initiation, propagation and junction in statics and dynamics. Computational Mechanics, 2007, 40(3): 473–495

[58] Zi G, Rabczuk T, Wall W. Extended meshfree methods without branch enrichment for cohesive cracks. Computational Mechanics, 2007, 40(2): 367–382

[59] Rabczuk T, Zi G. A meshfree method based on the local partition of unity for cohesive cracks. Computational Mechanics, 2007, 39(6): 743–760

[60] Rabczuk T, Areias P. A meshfree thin shell for arbitrary evolving cracks based on an extrinsic basis. Computer Modeling in Engineering & Sciences, 2006, 16(2): 115–130

[61] Rabczuk T, Zi G, Bordas S, Nguyen-Xuan H. A geometrically non-linear three-dimensional cohesive crack method for reinforced concrete structures. Engineering Fracture Mechanics, 2008, 75(16): 4740–4758

[62] Rabczuk T, Gracie R, Song J H, Belytschko T. Immersed particle method for fluid-structure interaction. International Journal for Numerical Methods in Engineering, 2010, 81(1): 48–71

[63] Rabczuk T, Bordas S, Zi G. On three-dimensional modelling of crack growth using partition of unity methods. Computers & Structures, 2010, 88(23–24): 1391–1411

[64] Amiri F, Anitescu C, Arroyo M, Bordas S, Rabczuk T. XLME interpolants, a seamless bridge between XFEM and enriched meshless methods. Computational Mechanics, 2014, 53(1): 45–57

[65] Amiri F, Millán D, Shen Y, Rabczuk T, Arroyo M. Phase-field modeling of fracture in linear thin shells. Theoretical and Applied Fracture Mechanics, 2014, 69: 102–109

[66] Hughes T J, Cottrell J A, Bazilevs Y. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 2005, 194(39–41): 4135–4195

[67] Cottrell J A, Hughes T J, Bazilevs Y. Isogeometric Analysis: Toward Integration of CAD and FEA. New York: John Wiley & Sons, 2009

[68] Benson D, Bazilevs Y, Hsu M C, Hughes T. Isogeometric shell analysis: The Reissner-Mindlin shell. Computer Methods in Applied Mechanics and Engineering, 2010, 199(5–8): 276–289

[69] Benson D, Bazilevs Y, Hsu M C, Hughes T. A large deformation, rotation-free, isogeometric shell. Computer Methods in Applied Mechanics and Engineering, 2011, 200(13–16): 1367–1378

[70] Nguyen-Thanh N, Zhou K, Zhuang X, Areias P, Nguyen-Xuan H, Bazilevs Y, Rabczuk T. Isogeometric analysis of large-deformation thin shells using RHT-splines for multiple-patch coupling. Computer Methods in Applied Mechanics and Engineering, 2017, 316: 1157–1178

[71] Nguyen-Thanh N, Kiendl J, Nguyen-Xuan H, Wüchner R, Bletzinger K U, Bazilevs Y, Rabczuk T. Rotation free isogeometric thin shell analysis using PHT-splines. Computer Methods in Applied Mechanics and Engineering, 2011, 200(47–48): 3410–3424

[72] Pawar A, Zhang Y, Jia Y, Wei X, Rabczuk T, Chan C, Anitescu C. Adaptive FEM-based nonrigid image registration using truncated hierarchical B-splines. Computers & Mathematics with Applications (Oxford, England), 2016, 72(8): 2028–2040

[73] Jia Y, Zhang Y, Rabczuk T. A novel dynamic multilevel technique for image registration. Computers & Mathematics with Applications (Oxford, England), 2015, 69(9): 909–925

[74] Thai T, Rabczuk T, Bazilevs Y, Meschke G. A higher-order stress-based gradient-enhanced damage model based on isogeometric analysis. Computer Methods in Applied Mechanics and Engineering, 2016, 304: 584–604

[75] Ghasemi H, Park H S, Rabczuk T. A multi-material level set-based topology optimization of flexoelectric composites. Computer Methods in Applied Mechanics and Engineering, 2018, 332: 47–62

[76] Ghasemi H, Park H, Rabczuk T. A level-set based IGA formulation for topology optimization of flexoelectric materials. Computer Methods in Applied Mechanics and Engineering, 2017, 313: 239–258

[77] Ghasemi H, Park H, Rabczuk T. Multi-material level set-based topology optimization of flexoelectric composites. Computer Methods in Applied Mechanics and Engineering, 2018, 332: 47–62

[78] Ghasemi H, Brighenti R, Zhuang X, Muthu J, Rabczuk T. Optimal fiber content and distribution in fiber-reinforced solids using a reliability and NURBS based sequential optimization approach. Structural and Multidisciplinary Optimization, 2015, 51(1): 99–112

[79] Anitescu C, Hossain M, Rabczuk T. Recovery-based error estimation and adaptivity using high-order splines over hierarchical T-meshes. Computer Methods in Applied Mechanics and Engineering, 2018, 328: 638–662

[80] Nguyen-Thanh N, Nguyen-Xuan H, Bordas S, Rabczuk T. Isogeometric analysis using polynomial splines over hierarchical T-meshes for two-dimensional elastic solids. Computer Methods in Applied Mechanics and Engineering, 2011, 200(21–22): 1892–1908

[81] Nguyen V, Anitescu C, Bordas S, Rabczuk T. Isogeometric analysis: An overview and computer implementation aspects. Mathematics and Computers in Simulation, 2015, 117: 89–116

[82] Anitescu C, Jia Y, Zhang Y, Rabczuk T. An isogeometric collocation method using super-convergent points. Computer Methods in Applied Mechanics and Engineering, 2015, 284: 1073–1097

[83] Nguyen B, Tran H, Anitescu C, Zhuang X, Rabczuk T. An isogeometric symmetric Galerkin boundary element method for two-dimensional crack problems. Computer Methods in Applied Mechanics and Engineering, 2016, 306: 252–275

[84] Jia Y, Anitescu C, Ghorashi S, Rabczuk T. Extended isogeometric analysis for material interface problems. IMA Journal of Applied Mathematics, 2015, 80(3): 608–633

[85] Ghorashi S, Valizadeh N, Mohammadi S, Rabczuk T. T-spline based XIGA for fracture analysis of orthotropic media. Computers & Structures, 2015, 147: 138–146

[86] Nguyen-Thanh N, Valizadeh N, Nguyen M, Nguyen-Xuan H, Zhuang X, Areias P, Zi G, Bazilevs Y, De Lorenzis L, Rabczuk T. An extended isogeometric thin shell analysis based on Kirchhoff-Love theory. Computer Methods in Applied Mechanics and Engineering, 2015, 284: 265–291

[87] Chan C, Anitescu C, Rabczuk T. Volumetric parametrization from a level set boundary representation with PHT-splines. Computer Aided Design, 2017, 82: 29–41

[88] Pérez N, Carbonari R C, Andrade M A B, Buiochi F, Adamowski J C. A FEM-based method to determine the complex material properties of piezoelectric disks. Ultrasonics, 2014, 54(6): 1631–1641

[89] Kiyono C Y, Pérez N, Silva E C N. Determination of full piezoelectric complex parameters using gradient-based optimization algorithm. Smart Materials and Structures, 2016, 25(2): 025019

[90] Pérez N, Buiochi F, Brizzotti Andrade M, Adamowski J. Numerical characterization of piezoceramics using resonance curves. Materials (Basel), 2016, 9(2): 71

[91] Vu-Bac N, Duong T, Lahmer T, Zhuang X, Sauer R, Park H, Rabczuk T. A NURBS-based inverse analysis for reconstruction of nonlinear deformations in thin shell structures. Computer Methods in Applied Mechanics and Engineering, 2018, 331: 427–455

[92] Ikeda T. Piezoelectricity. Oxford: Oxford University Press, 1990

[93] Yang J. An introduction to the theory of piezoelectricity. Vol. 9. Berlin: Springer, 2004

[94] Allik H, Hughes T J. Finite element method for piezoelectric vibration. International Journal for Numerical Methods in Engineering, 1970, 2(2): 151–157

[95] Holland R. Representation of dielectric, elastic, and piezoelectric losses by complex coefficients. IEEE Transactions on Sonics and Ultrasonics, 1967, 14(1): 18–20

[96] Lahrner T, Kaltenbacher M, Kaltenbacher B, Lerch R, Leder E. FEM-based determination of real and complex elastic, dielectric, and piezoelectric moduli in piezoceramic materials. IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 2008, 55(2): 465–475

[97] Huyer W, Neumaier A. Global optimization by multilevel coordinate search. Journal of Global Optimization, 1999, 14(4): 331–355

[98] Vu-Bac N, Lahmer T, Zhuang X, Nguyen-Thoi T, Rabczuk T. A software framework for probabilistic sensitivity analysis for computationally expensive models. Advances in Engineering Software, 2016, 100: 19–31

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature