Abstract
Optimal design of piezoelectric vibration devices used in medical fields is a typical inequality constrained optimization problem. Constrained variable method, as the best algorithm to solve nonlinear constraint optimization problem is used in the optimal design of the piezoelectric devices. The object functions and constrain functions are automatically calculated with Finite Element Analysis software ANSYS, and the results will return to the main optimization program. Results show that compared with the given initial parameters, the vibration amplitude of the optimum piezoelectric vibration device is dramatically increased under the constraint that the stress concentration of the cutting tool is within limit. Experiments are also performed to measure the vibration amplitude. Results are in great agreement with the theoretical one.
1. Introduction
Piezoelectric vibration devices are widely used in cleaning, machining, dental and surgery. When apply in medical field, they must meet the requirements of high output power, high efficiency, lightweight and high reliability [1]. For piezoelectric devices used in bone surgery, the demand of high power and high reliability are incompatible. For example, change of structure is the necessary method to achieve high output power, however this will also lead to stress concentration and increase the probability of device breakage [23]. In a word, design of piezoelectric vibration devices is an optimization problem with multidesign variants, multiconstraints and multiobject.
Previous studies of piezoelectric vibration devices optimization are mostly based on electrical equivalent circuit or based on experience of the designers [46]. Lockwood etc. [7] describes a method of modelling transducers using network theory, which is determined from measurement of the transducer impedance in water and the pulseecho response of the system for a given electrical source and load. P. Tierce etc. [8] pointed out that finiteelement method is a new way to predict vibratory behaviour and can provide information for optimization of piezoelectric transducers. Pascal [9] proposed a method which integrated finiteelement method and some experimental results to study the effects of elementary ceramic rods widthtothickness ratio.
Among numerous optimization algorithms, Constrained Variable Method is an excellent nonlinear constrained optimization method, which was firstly introduced by Wilson in 1963 and consummated by Han in 1977 [10] and Powell in 1978 [11]. This method has the advantage of fast convergence speed, high efficiency, good stability and excellent adaptability. Thus, it has been widely used in engineering applications [12]. Finite Element Analysis (FEA) is widely used in engineering field [13] and can be used to achieve some constraint variable.
In this paper, a new method to optimize piezoelectric vibration devices is proposed. This method integrates Constrained Variable Algorithm with (FEA) method, in which FEA will provide the value of optimization variant, constraint function and object function to the main algorithm program.
The remainder of the paper is organized as follows. Section 2 introduced the theory of Constrained Variable Method. Section 3 describes the proposed optimization method with its flow diagram. Section 4 presents the piezoelectric coupled FEA theory. Section 5 uses an example to illustrate the optimization procedures, and the effectiveness of the proposed methodology.
2. Constrained variable optimization method
An inequality constrained optimization problem can be described as:
$\mathrm{m}\mathrm{i}\mathrm{n}f\left(x\right),$
$\text{s.t.}{g}_{j}\left(x\right)\ge 0,j=1,\mathrm{}2,\dots ,m,$
${\underset{\_}{x}}_{i}\le {x}_{i}\le {\stackrel{}{x}}_{i},i=1,\mathrm{}2,\dots ,n,$
where $x$ is design variable vector, $g\left(x\right)$ is constrain function and is object function. ${\underset{\_}{x}}_{i}$ and ${\stackrel{}{x}}_{i}$ are the minimal and maximum limit of design variables.
To solve this optimization problem by constrained variable method, Eq. (1) should firstly convert to its quadratic programming subproblem, and then looks for the optimal results by Newton iteration method.
The corresponding quadratic programming problem of Eq. (1) is:
$\mathrm{m}\mathrm{i}\mathrm{n}\frac{1}{2}{d}^{T}{B}^{k}d+\nabla f{\left({x}^{k}\right)}^{T}d,$
$\text{s.t.}{g}_{j}\left({x}^{k}\right)+\nabla {g}_{j}{\left({x}^{k}\right)}^{T}d\ge 0,j=1,\mathrm{}2,\dots ,m,$
${\underset{\_}{x}}_{i}\le {x}_{i}^{k}+d\le {\stackrel{}{x}}_{i},i=1,\mathrm{}2,\dots ,n,$
where $d$ is search direction, $\nabla f\left({x}^{k}\right)$ and $\nabla {g}_{j}\left({x}^{k}\right)$ are gradient of object function and constraint functions, ${B}^{k}$ is revision Hessian Matrix of Lagrange function.
To solve Eq. (2) by Newton method, search direction must be firstly determined. Then search step can be calculated with onedimensional search method.
Zanggwill penalty function is used as onedimensional search function, which is:
$+\sum _{i=1}^{n}\lambda {\text{\hspace{0.17em}}}_{i}\mathrm{m}\mathrm{a}\mathrm{x}\left\{\text{\hspace{0.17em}}0,{x}_{i}^{k}+\alpha \text{\hspace{0.17em}}{d}^{k}{\stackrel{}{x}}_{i},{\underset{\_}{x}}_{i}{x}_{i}^{k}\alpha \text{\hspace{0.17em}}{d}^{k}\right\},$
where $\alpha $ is search step, to increase convergence speed, imprecise search is usually adopted, the searching rule is:
where $\mathrm{\Phi}\left(\alpha \right)$ is a segment function. Firstly let $\alpha =\text{1}$, then calculate Eq. (4), if it is not satisfied, let $\alpha $ equal to 0.1, 0.01, … until Eq. (4) is satisfied.
To solve Eq. (1), the derivatives of object function and constraint function to the design variables should be calculated, this process is called characteristic sensitivity analysis. Finite difference method is adopted here to calculate the sensitivity:
where $u$ is object function or constraint function, $\mathrm{\Delta}{x}_{j}$ is step size, and can be selected by experience.
3. The proposed optimization method
Flow diagram of the proposed optimization method is shown in Fig. 1.
Fig. 1Flow diagram of the proposed optimization method
Step 1: Give an initial value ${x}^{0}$ and convergence accuracy ${\epsilon}_{1}$, ${\epsilon}_{2}$, let ${B}^{0}=I$, $k=\text{0}$.
Step 2: Call interface subprogram PS1, calculate the constraint function $g\left({x}^{k}\right)$ and object function $f\left({x}^{k}\right)$ with FEA, and also the gradient of constraint function $\nabla g\left({x}^{k}\right)$ and object function $\nabla f\left({x}^{k}\right)$.
Step 3: Form quadratic programming problem, obtain search direction ${d}^{k}$.
Step 4: Perform onedimensional search, and get ${\alpha}^{k}\text{,}$${x}^{k+1}={x}^{k}+{\alpha}^{k}{d}^{k}\text{,}$ call interface subprogram PS1 to calculate $g\left({x}^{k+1}\right)$, $f\left({x}^{k+1}\right)$, $\nabla g\left({x}^{k+1}\right)$ and $\nabla f\left({x}^{k+1}\right)$.
Step 5: Judge whether the convergence criteria is satisfied. If the answer is yes, optimization result is obtained, otherwise go to step 6.
Step 6: Revise variable metric matrix ${B}_{k}$. Let $k=k+1$ and go to step 3.
The function of interface subprogram PS1 may transmits data between optimization program and ANSYS software via data file called dat1. These data include design variables, constrain variables and object variables.
Dat1 is a batch file programed by ANSYS APDL. The main function of this file includes: (1) Read the design variables, constrain variables and object variables from main program. (2) Parametric modeling and mesh automatically. (3) Loading and solve the problem. (4) Post processing and extract the needed data. (5) Output the constrain variables and object variables to the main program.
4. Piezoelectric coupled FEA theory
A typical structure of piezoelectric vibration devices is shown in Fig. 2, which mainly includes the bolt, electrode, piezoelectric ceramics, horn and cutting tool. Based on piezoelectric effect, the ceramic convert electric signal into mechanical vibration at the tip of the cutting tool, which vibrates longitudinally to scrape or cut bone.
Fig. 2Typical structure of piezoelectric vibration devices
Assumptions are made in Finite Element Analysis, which are: (1) Thickness of bonding agent between piezoelectric ceramic plates is negligible. (2) Thickness of the electrode slice is negligible compared with piezoelectric ceramic plates. (3) Mechanical friction between piezoelectric ceramic and other parts is negligible. (4) Performance loss of piezoelectric material and other parts is negligible. (5) Piezoelectric material is isotropic, polarization direction is $Z$ axis.
Fundamental equations of piezoelectric vibration devices include constitutive equation, motion equation, geometric equation and boundary conditions.
4.1. Constitutive equation
4.1.1. Piezoelectric constitutive equation
Suppose that the piezoelectric polarization direction, electric field and mechanical vibration direction are all parallel to $Z$ axis. The piezoelectric constitutive equation is:
where $T$, $S$, $D$ and $E$ are stress vector, strain vector, electrical displacement vector and electrical field intensity vector. ${c}_{33}^{E}$ is elastic stiffness coefficient under constant electrical field intensity, ${e}_{33}$ is piezoelectric stress constant and ${\epsilon}_{33}^{S}$ is material dielectric constant.
4.1.2. Constitutive equation for other parts
Materials of other part of the piezoelectric device are all isotropic linear elastic materials, their constitutive equation is:
where $\lambda $ and $G$ are constant parameters:
where $E$ and $\nu $ are Young modulus and Poisson’s ratio.
4.2. Motion equation
The basic equation of motion of elastimechanics under cylindrical coordinate is:
where ${F}_{r}$ and ${F}_{z}$ are body forces of axial direction and $Z$ direction. ${u}_{r}$ and $w$ are displacement component of axial direction and $Z$ direction.
4.3. Geometric equation
The relationship between stain and displacement of elastic material, which is also called geometric equation, is given under cylindrical coordinate:
4.4. Boundary conditions
Boundary at 1 and 8 location (Fig. 2) is freesurface condition; the stress at the surface is zero, i.e. ${\sigma}_{ij}{n}_{i}=\text{0}$. Boundary of 2, 4 and 6 are electrical free surface, the electric density at electrical charge free surface is zero, which is ${D}_{i}{n}_{j}=\text{0}$. The electric potential at electrical potential free surface is zero, which is $\varphi =\text{0}$. At boundary of 3 and 5, given electrical charge free surface the electric charge density boundary condition is ${D}_{i}{n}_{j}=\stackrel{}{w}$. Given electrical potential free surface the electric potential boundary condition is $\varphi =\stackrel{}{\varphi}$. Location 7 is clamped boundary and ${u}_{i}=\text{0}$, $i=x$, $y$, $z$, where ${u}_{i}$ is the displacement at three axis.
4.5. Finite element method
Constitutive Eq. (6)(8), motion Eq. (19), geometric Eq. (10) and boundary conditions constitute the governing equation of the piezoelectric device and can be solved by numerical method, such as finite element method.
With shape function, displacement and electric potential difference of an element is:
where $\text{{}{u}_{c}\text{}}$ is element displacement vector, ${V}_{c}$ is element voltage vector, $\text{{}u\text{}}$ is node displacement matrix, $\text{{}V\text{}}$ is node electromotive force vector, $\left[{N}^{u}\right]$ is displacement shape function matrix, $\left[{N}^{V}\right]$ is voltage shape function matrix:
where ${N}_{i}$ is shape function of node $i$:
Motion equation can convert into linear form:
where $\left[{M}_{e}\right]$ is mass matrix:
where ${V}_{e}$ is the volume of the elastic body, which includes mental material and piezoelectric material, $\rho $ is material density, and $\left[{K}_{e}\right]$ is stiffness matrix:
where $\left[C\right]$ is elastic constant matrix of materials:
where $\left[{C}_{e}\right]$ is damping matrix:
where $\alpha $ and $\beta $ are Rayleigh damping constant, $\left[{K}_{e}^{d}\right]$ is dielectric matrix and $\left[{K}_{e}^{c}\right]$ is piezoelectric coupling matrix:
where ${V}_{p}$ is the volume of piezoelectric material, $\left[\mu \right]$ is dielectric constant matrix of piezoelectric material, and:
where $\left[e\right]$ is piezoelectric matrix.
5. A case study
5.1. Piezoelectric device and FEA modeling
The piezoelectric device in this case study is used in bone surgery, which is shown in Fig. 3.
Fig. 3Structure of the piezoelectric device
Some of the structure sizes are determined by analytic or empirical calculation, as shown in Table 1.
Table 1Parameters of the structure
Parameters  ${L}_{1}$  ${L}_{2}$  ${L}_{3}$  ${L}_{5}$  ${L}_{6}$  ${L}_{8}$  ${D}_{1}$  ${D}_{3}$ 
Value (mm)  12  12  15  12  9  12  15  ${D}_{2}$+3 
Other parameters, include ${L}_{4}$, ${L}_{7}$, ${D}_{2}$, ${D}_{4}$ are design variables. The optimal object is large vibration amplitude at the tip of the cutting tool, and constrain function is the stress at transitional part and the teeth root (shown in Fig. 3) should not exceed their empirical limit, in order to avoid the cutting tool crack due to stress concentration. Thus, the mathematics model of the piezoelectric device optimization problem is given as:
$\mathrm{min}{f}_{1}\left(x\right),$
$\text{s.t.}$
${g}_{1}\left(x\right)=1\frac{{\sigma}_{eqv}^{A}}{\left[{\sigma}_{eqv}^{}\right]}\ge 0,$
${g}_{2}\left(x\right)=1\frac{{\sigma}_{eqv}^{B}}{\left[{\sigma}_{eqv}^{}\right]}\ge 0,$
$5\le {L}_{4}\le 8,$
$3\le {L}_{7}\le 6,$
$8\le {D}_{2}\le 15,$
$5\le {D}_{4}\le 12,$
where ${f}_{1}\left(x\right)$ is vibration amplitude at the tip of the cutting tool, ${g}_{1}\left(x\right)$ and ${g}_{2}\left(x\right)$ are the constrain functions to avoid crack at transitional part and tooth root. ${\sigma}_{eqv}^{A}$ and ${\sigma}_{eqv}^{B}$ are the equivalent stress at transitional part and teeth root, $\left[{\sigma}_{eqv}^{}\right]$ is the empirical stress limit at transitional part and teeth root. ${f}_{1}\left(x\right)$, ${\sigma}_{eqv}^{A}$ and ${\sigma}_{eqv}^{B}$ are all calculated by commercial FEA software ANSYS.
In ANSYS, properties of steel, titanium, and PZT8 ceramic are listed in Table 2.
Table 2Material properties
Materials  Young’s modulus (N/m)  Density (kg/m^{3})  Poisson’s ratio 
Steel  2.1×1011  7896  0.29 
Titanium  1.05×1011  4450  0.33 
PZT8  –  7500  – 
Piezoelectric rings are modelled using 3D coupledfield solid elements SOLID5. The dielectric matrix, piezoelectric matrix (unit: C/m^{2}) and elastic constant matrix (unit: 1010 N/m^{2}) of PZT8 are:
$\left[e\right]=\left[\begin{array}{ccc}0& 0& 4.7\\ 0& 0& 4.7\\ 0& 0& 12.9\\ 0& 10.3& 0\\ 10.3& 0& 0\\ 0& 0& 0\end{array}\right],$
$\left[c\right]=\left[\begin{array}{cccccc}15.6& 8.86& 8.45& 0& 0& 0\\ & 15.6& 8.45& 0& 0& 0\\ & & 13.0& 0& 0& 0\\ & & & 3.13& 0& 0\\ & & & & 3.13& 0\\ & & & & & 3.38\end{array}\right].$
Modal analysis is firstly performed using the Lancoze method. Then harmonic analysis is performed using Sparse Matrix Solver. The input voltage varied sinusoidally between +/–600, which is shown in Fig. 4.
Fig. 4Loads apply on piezoelectric ceramic rings
5.2. Results and discussion
The initial value of design variables ${L}_{4}$, ${L}_{7}$, ${D}_{2}$ and ${D}_{4}$ are 6 mm, 4 mm, 12 mm and 8 mm. After five iterations and call the FEA software 31 times, the optimal results are achieved, which is ${L}_{4}=\text{6.78}\text{}\text{mm}$, ${L}_{7}=\text{3.54}\text{}\text{mm}$, ${D}_{2}=\text{10.57}\text{}\text{mm}$, ${D}_{4}=\text{6.26}\text{}\text{mm}$.
The values of design variables, constraint functions and object function of each iteration are shown in Table 3.
Table 3Optimization process variables
Parameters  Iteration  
0  1  2  3  4  5  
${L}_{4}\text{(mm)}$  6.0000  6.2353  6.4172  6.5356  6.7791  6.7833 
${L}_{7}\text{(mm)}$  4.0000  4.0255  3.9828  3.6271  3.5214  3.5421 
${D}_{2}\text{(mm)}$  12.0000  12.1327  11.8201  11.1529  10.6282  10.5719 
${D}_{4}\text{(mm)}$  8.0000  8.0491  7.8934  6.9453  6.3185  6.2624 
${g}_{1}\left(x\right)$  0.1749  0.1828  0.1777  0.1049  0.0732  0.0805 
${g}_{2}\left(x\right)$  0.4371  0.4342  0.4318  0. 4263  0.4256  0.4264 
${f}_{1}\left(x\right)\text{}\text{(mm)}$  0.0973  0.0939  0.1127  0.1471  0.1482  0.1514 
From Table 3, after five iterations, constraint function ${g}_{1}\left(x\right)$ is close to zero, that means the design variables cannot adjust otherwise ${g}_{1}\left(x\right)$ will be negative and violates the constraint condition, and that point is the optimal results. At the beginning, ${f}_{1}\left(x\right)$ is 97.3 um, after the iteration ${f}_{1}\left(x\right)$ is 151.4 um, the optimization object is achieved.
After optimization, equivalent stress distribution of the cutting tool is relatively poorer compared with that of initial value. With initial design variant, the maximum equivalent stress at the cutting tool is 17.7 MPa, while after optimization, this value is 19.0 MPa.
Fig. 5Equivalent stress distribution before and after optimization: a) before optimization, b) after optimization
a)
b)
Piezoelectric vibration devices with optimum design parameters is manufactured and assembled, which is shown in Fig. 6.
Admittance measurement is carried out with dynamic signal analysis instrument SR785. Results in Fig. 7 show that the resonance frequency of this device is 40.382 kHz. The value is 40.350 kHz in FEA analysis.
Fig. 6Piezoelectric vibration devices with optimum design parameters
Fig. 7Measurement of resonance frequency of the optimal piezoelectric device
Vibration amplitude of the piezoelectric device is measured with two methods. One is highspeed imaging method with Keyence VHX100 digital microscope. Fig. 8 shows the schematic diagram of this measurement scheme. The measurement system is composed of a microscope with highspeed digital camera lens, which can pick up the image when piezoelectric device vibrates with ultrasonic frequency.
Fig. 9 shows the images taking with highspeed microscope, where one grid of the scale represents 100 µm. In Fig. 6(b), the vibration belt shows that displacement of the tip is one and a half grids, which is nearly 150 µm.
Vibration amplitude measured with highspeed imaging instrument is an approximate value, precise measurement can be achieved by laser Doppler vibrometer, which is shown in Fig. 10. This measurement system is composed of laser Doppler vibrometer (CV700), oscilloscope, and piezoelectric signal generator.
Fig. 8Measurement of tip displacement with highspeed imaging instrument
Fig. 9Image of the cutting tool of piezoelectric device: a) resting and b) vibrating condition
a)
b)
Output signal in Fig. 10 is vibration velocity, which will be integrated to obtain vibration displacement. Fig. 11 shows the vibration displacement at the tip of the piezoelectric device cutting tool, from which the vibration amplitude can be precisely measured, and the results is also near 150 um, consistent with the FEA method and imaging method.
Fig. 10Measurement of tip displacement with laser Doppler vibrometer
The optimized piezoelectric tool is used to cut different kinds of bones, which is shown in Fig. 12. Spongy bone is from the vertebra and compact bone is from the femur of a pig. The speed when cutting spongy bone is nearly 0.01 m/s, while the speed when cutting compact bone is slightly slower than cutting spongy bone. More than 30 clinical cases also show that the cutting speed of the optimized piezoelectric tool is higher than tradition electric bone grinding tool.
Fig. 11Vibration displacement of cutting tool tip
Fig. 12Cutting different bone tissues with optimized piezoelectric tool: a) spongy bone, b) compact bone
a)
b)
Lifetime of 85 piezoelectric bone tools that used in clinical cases are recorded, among these tools, 60 of them are optimized with constrained variable algorithm and FEA method, others are designed with traditional method or by experience of the designer. It is shown in Table 4, lifetime of the optimized tools are longer than those tools that not been optimized.
Table 4Lifetime of piezoelectric bone with and without design optimization
Lifetime (h)  $\text{1}<T\le \text{30}$  $\text{30}<T\le \text{50}$  $\text{50}<T\le \text{100}$  $T>\text{100}$  Total  Average lifetime (h) 
Tools with design optimization  1  8  29  22  60  ≈80 
Tools without design optimization  1  9  10  5  25  ≈65 
From the above analysis and experiments, the optimized piezoelectric tool can output higher energy, and more efficient in cutting bone, which means doctors will save labour and time to finish the surgery. Although local stress on the tool will increase, lifetime of piezoelectric tools is somewhat improved due to the shape optimization. Conclusions can be drawn that the optimization method with Constrained Variable Algorithm and FEA is very effective to design a high performance piezoelectric surgery tool.
6. Conclusions
Constrained Variable Algorithm and FEA method are used for the optimization of piezoelectric vibration device. Results show that for nonlinear optimization problem with multidesign variants and multiconstraints, Constrained Variable Algorithm has fast convergence speed. And when this algorithm used in engineering fields, FEA software is a good supplement for calculating constraint function and object function.
References

Dong Sun, Zhaoying Zhou, Liu Y. H., Shen Z. Development and application of ultrasonic surgical instruments. Transactions on Biomedical Engineering, Vol. 44, Issue 6, 1997, p. 462467.

Ying Chen, Zhou Z. Y. Effects of different tissue loads on high power ultrasonic surgery scalpel. Ultrasound in Medicine and Biology, Vol. 32, Issue 3, 2006, p. 415420.

Ying Chen, Rui Kang Design for reliability of ultrasonic bone knife. Advanced Science Letters, Vol. 4, Issue 67, 2011, p. 21472151.

Gu Yujiong, Zhou Zhaoying, Yao Jian Ultrasonic vibration system four terminal network method and its application. Chinese Journal of Mechanical Engineering, Vol. 33, Issue 3, 1997, p. 94101.

Sergio S. Design of an optimized highpower ultrasonic transducer. Muhlen Ultrasonics Symposium, 1990, p. 16311633.

Eisner E., Seager J. S. A longitudinally resonant stub for vibrations of large amplitude. Ultrasonics, 1965, p. 8898.

Lockwood Geoffrey R., Foster F. Stuart Modeling and optimization of highfrequency ultrasound transducers. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, Vol. 41, Issue 2, 1994, p. 225230.

Tierce P., Decarpigny J. N., Hamonic B., Decarpigny J. N. Power Sonic and Ultrasonic Transducer Design. Springer, Heidelberg, 1988.

Pascal C. Optimizing ultrasonic transducers based on piezoelectric composites using a finiteelement method. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, Vol. 37, Issue 2, 1990, p. 135140.

Han S. P. A globally convergent method for nonlinear programming. Journal of Optimization Theory and Applications, Vol. 22, Issue 3, 1977, p. 297309.

Powell M. J. D. Algorithms for nolinear constraints that use lagrangian function. Mathematical Programming, Vol. 14, 1980, p. 224248.

Bartholomew M. Nonlinear optimization with engineering applications. Springer Science Business Media, LLC, 2008.

Ying Chen, Bingdong Liu, Rui Kang Study of vibratory stress relief effect of quartz flexible accelerometer with FEA method. Journal of Vibroengineering, Vol. 15, Issue 2, 2013, p. 784793.