Published: 30 September 2014

Optimal design of piezoelectric vibration devices with constrained variable algorithm and FEA

Ying Chen1
Cui Ye2
Rui Kang3
1, 2, 3Science and Technology on Reliability and Environmental Engineering Laboratory, Beihang University, Beijing, 100191, P. R. China
Corresponding Author:
Ying Chen
Views 58
Reads 25
Downloads 1084


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 [2-3]. In a word, design of piezoelectric vibration devices is an optimization problem with multi-design variants, multi-constraints and multi-object.

Previous studies of piezoelectric vibration devices optimization are mostly based on electrical equivalent circuit or based on experience of the designers [4-6]. 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 pulse-echo response of the system for a given electrical source and load. P. Tierce etc. [8] pointed out that finite-element 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 finite-element method and some experimental results to study the effects of elementary ceramic rods width-to-thickness ratio.

Among numerous optimization algorithms, Constrained Variable Method is an excellent non-linear 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:

s.t. gjx0, j=1, 2,,m,
x_ixix-i, i=1, 2,,n,

where x is design variable vector, g(x) is constrain function and is object function. x_i and 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 sub-problem, and then looks for the optimal results by Newton iteration method.

The corresponding quadratic programming problem of Eq. (1) is:

s.t. gjxk+gjxkTd0, j=1, 2,,m,
x_ixik+dx-i, i=1, 2,,n,

where d is search direction, f(xk) and gj(xk) are gradient of object function and constraint functions, Bk 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 one-dimensional search method.

Zanggwill penalty function is used as one-dimensional search function, which is:


where α is search step, to increase convergence speed, imprecise search is usually adopted, the searching rule is:


where Φ(α) is a segment function. Firstly let α=1, then calculate Eq. (4), if it is not satisfied, let α 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:

uixj=ui(x+Δxj)-ui(x)Δxj, j=1, 2,,n,

where u is object function or constraint function, Δxj 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

Flow diagram of the proposed optimization method

Step 1: Give an initial value x0 and convergence accuracy ε1, ε2, let B0=I, k=0.

Step 2: Call interface subprogram PS1, calculate the constraint function g(xk) and object function f(xk) with FEA, and also the gradient of constraint function g(xk) and object function f(xk).

Step 3: Form quadratic programming problem, obtain search direction dk.

Step 4: Perform one-dimensional search, and get αk,xk+1=xk+αkdk, call interface subprogram PS1 to calculate g(xk+1), f(xk+1), g(xk+1) and f(xk+1).

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 Bk. 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

Typical 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. c33E is elastic stiffness coefficient under constant electrical field intensity, e33 is piezoelectric stress constant and ε33S 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 λ and G are constant parameters:

λ=Eν(1+ν)(1-2ν), G=E2(1+ν),

where E and ν are Young modulus and Poisson’s ratio.

4.2. Motion equation

The basic equation of motion of elasti-mechanics under cylindrical coordinate is:


where Fr and Fz are body forces of axial direction and Z direction. ur 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:

εr=urr, εθ=urr, εz=wz, γrz=urz+wr.

4.4. Boundary conditions

Boundary at 1 and 8 location (Fig. 2) is free-surface condition; the stress at the surface is zero, i.e. σijni=0. Boundary of 2, 4 and 6 are electrical free surface, the electric density at electrical charge free surface is zero, which is Dinj=0. The electric potential at electrical potential free surface is zero, which is ϕ=0. At boundary of 3 and 5, given electrical charge free surface the electric charge density boundary condition is Dinj=w-. Given electrical potential free surface the electric potential boundary condition is ϕ=ϕ-. Location 7 is clamped boundary and ui=0, i=x, y, z, where ui 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:

uc=NuTu, Vc=NVTV,

where {uc} is element displacement vector, Vc is element voltage vector, {u} is node displacement matrix, {V} is node electromotive force vector, [Nu] is displacement shape function matrix, [NV] is voltage shape function matrix:


where Ni is shape function of node i:


Motion equation can convert into linear form:


where [Me] is mass matrix:


where Ve is the volume of the elastic body, which includes mental material and piezoelectric material, ρ is material density, and [Ke] is stiffness matrix:


where C is elastic constant matrix of materials:


where [Ce] is damping matrix:


where α and β are Rayleigh damping constant, [Ked] is dielectric matrix and [Kec] is piezoelectric coupling matrix:


where Vp is the volume of piezoelectric material, [μ] is dielectric constant matrix of piezoelectric material, and:


where [e] 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

Structure 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

Value (mm)

Other parameters, include L4, L7, D2, D4 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:

g1(x)=1-σeqvAσeqv 0,
g2x=1-σeqvBσeqv 0,

where f1(x) is vibration amplitude at the tip of the cutting tool, g1(x) and g2(x) are the constrain functions to avoid crack at transitional part and tooth root. σeqvA and σeqvB are the equivalent stress at transitional part and teeth root, [σeqv ] is the empirical stress limit at transitional part and teeth root. f1(x), σeqvA and σeqvB 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

Young’s modulus (N/m)
Density (kg/m3)
Poisson’s ratio

Piezoelectric rings are modelled using 3-D coupled-field solid elements SOLID5. The dielectric matrix, piezoelectric matrix (unit: C/m2) and elastic constant matrix (unit: 1010 N/m2) of PZT8 are:


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

Loads apply on piezoelectric ceramic rings

5.2. Results and discussion

The initial value of design variables L4, L7, D2 and D4 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 L4=6.78 mm, L7=3.54 mm, D2=10.57 mm, D4=6.26 mm.

The values of design variables, constraint functions and object function of each iteration are shown in Table 3.

Table 3Optimization process variables

L4 (mm)
L7 (mm)
D2 (mm)
D4 (mm)
0. 4263
f1x (mm)

From Table 3, after five iterations, constraint function g1(x) is close to zero, that means the design variables cannot adjust otherwise g1(x) will be negative and violates the constraint condition, and that point is the optimal results. At the beginning, f1(x) is 97.3 um, after the iteration f1(x) 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

Equivalent stress distribution before and after optimization:  a) before optimization, b) after optimization


Equivalent stress distribution before and after optimization:  a) before optimization, b) after optimization


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

Piezoelectric vibration devices  with optimum design parameters

Fig. 7Measurement of resonance frequency of the optimal piezoelectric device

Measurement of resonance frequency  of the optimal piezoelectric device

Vibration amplitude of the piezoelectric device is measured with two methods. One is high-speed imaging method with Keyence VHX-100 digital microscope. Fig. 8 shows the schematic diagram of this measurement scheme. The measurement system is composed of a microscope with high-speed digital camera lens, which can pick up the image when piezoelectric device vibrates with ultrasonic frequency.

Fig. 9 shows the images taking with high-speed 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 high-speed 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 high-speed imaging instrument

Measurement of tip displacement with high-speed imaging instrument

Fig. 9Image of the cutting tool of piezoelectric device: a) resting and b) vibrating condition

Image of the cutting tool of piezoelectric device: a) resting and b) vibrating condition


Image of the cutting tool of piezoelectric device: a) resting and b) vibrating condition


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

Measurement 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

Vibration displacement of cutting tool tip

Fig. 12Cutting different bone tissues with optimized piezoelectric tool: a) spongy bone, b) compact bone

Cutting different bone tissues with optimized piezoelectric tool: a) spongy bone, b) compact bone


Cutting different bone tissues with optimized piezoelectric tool: a) spongy bone, b) compact bone


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)
Average lifetime (h)
Tools with design optimization
Tools without design optimization

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 multi-design variants and multi-constraints, 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.


  • 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. 462-467.
  • 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. 415-420.
  • Ying Chen, Rui Kang Design for reliability of ultrasonic bone knife. Advanced Science Letters, Vol. 4, Issue 6-7, 2011, p. 2147-2151.
  • 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. 94-101.
  • Sergio S. Design of an optimized high-power ultrasonic transducer. Muhlen Ultrasonics Symposium, 1990, p. 1631-1633.
  • Eisner E., Seager J. S. A longitudinally resonant stub for vibrations of large amplitude. Ultrasonics, 1965, p. 88-98.
  • Lockwood Geoffrey R., Foster F. Stuart Modeling and optimization of high-frequency ultrasound transducers. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, Vol. 41, Issue 2, 1994, p. 225-230.
  • 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 finite-element method. IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control, Vol. 37, Issue 2, 1990, p. 135-140.
  • Han S. P. A globally convergent method for nonlinear programming. Journal of Optimization Theory and Applications, Vol. 22, Issue 3, 1977, p. 297-309.
  • Powell M. J. D. Algorithms for nolinear constraints that use lagrangian function. Mathematical Programming, Vol. 14, 1980, p. 224-248.
  • 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. 784-793.

About this article

12 April 2014
30 June 2014
30 September 2014
piezoelectric vibration device
design optimization
constrained variable algorithm