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

1Corresponding author

Journal of Vibroengineering, Vol. 16, Issue 6, 2014, p. 2962-2973.
Received 12 April 2014; received in revised form 4 June 2014; accepted 30 June 2014; published 30 September 2014

Copyright © 2014 JVE International Ltd. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Creative Commons License
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.

Keywords: piezoelectric vibration device, design optimization, FEA, constrained variable algorithm.

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:

(1)
x = x 1 , x 2 , , x n T ,
m i n f x ,
s.t.  g j x 0 ,       j = 1 ,   2 , , m ,
x _ i x i x - 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:

(2)
d = d 1 , d 2 , , d n T ,
m i n 1 2 d T B k d + f ( x k ) T d ,
s.t.  g j x k + g j x k T d 0 ,       j = 1 ,   2 , , m ,
x _ i x i k + d x - 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:

(3)
Φ α = f x k + α d k + j = 1 m μ j m a x 0 , - g j x k + α d k
          + i = 1 n λ i m a x 0 , x i k + α d k - x - i , x _ i - x i k - α d k ,

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

(4)
Φ α k Φ 0 + 0.1 α k Φ ' 0 ,

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:

(5)
u i x j = u i ( x + Δ x j ) - u i ( x ) Δ x j ,         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. 1. Flow 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. 2. Typical 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:

(6)
D = ε s s s E + e s s S , T = - e s s E + c s s E S ,

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:

(7)
σ r = λ θ + 2 G ε r , σ θ = λ θ + 2 G ε θ , σ z = λ θ + 2 G ε z , τ r z = G γ r z ,

where λ and G are constant parameters:

(8)
λ = E ν ( 1 + ν ) ( 1 - 2 ν ) ,       G = E 2 ( 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:

(9)
σ r r + τ r z z + σ r - σ θ r + F r = ρ 2 u r t 2 , τ r z r + σ z z + τ r z r + F z = ρ 2 w t 2 ,                    

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:

(10)
ε r = u r r ,       ε θ = u r r ,       ε z = w z ,       γ r z = u r z + w r .

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:

(11)
u c = N u T u ,       V c = N V T V ,

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:

(12)
[ N u ] T = N 1 0 0 N n 0 0 0 N 1 0 0 N n 0 0 0 N 1 0 0 N n ,

where Ni is shape function of node i:

(13)
[ N V ] T = N 1 N 2 N n .

Motion equation can convert into linear form:

(14)
M e 0 0 0 u ¨ V ¨ + C e 0 0 0 u ˙ V ˙ + K e K e c K e c T K e d u V = F e Q e ,

where [Me] is mass matrix:

(15)
M e = V e ρ N u N u T d V ,

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

(16)
K e = V e B u T C B u d V ,

where C is elastic constant matrix of materials:

(17)
B u = x 0 0 0 y 0 0 0 z x x 0 0 y z x 0 z N u T ,

where [Ce] is damping matrix:

(18)
C e = α M e + β K e ,

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

(19)
K e d = - V p B V T μ B V d V ,

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

(20)
B V = x y z N V T ,
(21)
K e c = V p B u T e B V d V ,

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. 3. Structure 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 1. Parameters 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 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:

x * = L 4 * , L 7 * , D 2 * , D 4 * T ,
min f 1 x ,
s.t.
g 1 ( x ) = 1 - σ e q v A σ e q v     0 ,
g 2 x = 1 - σ e q v B σ e q v   0 ,
5 L 4 8 ,
3 L 7 6 ,
8 D 2 15 ,
5 D 4 12 ,

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 2. Material properties

Materials
Young’s modulus (N/m)
Density (kg/m3)
Poisson’s ratio
Steel
2.1×1011
7896
0.29
Titanium
1.05×1011
4450
0.33
PZT8
7500

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:

μ = 7.97 0 0 0 7.97 0 0 0 5.13 × 10 - 9 ,
e = 0 0 - 4.7 0 0 - 4.7 0 0 12.9 0 10.3 0 10.3 0 0 0 0 0 ,
c = 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 .

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. 4. Loads 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 3. Optimization process variables

Parameters
Iteration
0
1
2
3
4
5
L 4   (mm)
6.0000
6.2353
6.4172
6.5356
6.7791
6.7833
L 7   (mm)
4.0000
4.0255
3.9828
3.6271
3.5214
3.5421
D 2   (mm)
12.0000
12.1327
11.8201
11.1529
10.6282
10.5719
D 4   (mm)
8.0000
8.0491
7.8934
6.9453
6.3185
6.2624
g 1 ( x )
0.1749
0.1828
0.1777
0.1049
0.0732
0.0805
g 2 ( x )
0.4371
0.4342
0.4318
0. 4263
0.4256
0.4264
f 1 x   (mm)
0.0973
0.0939
0.1127
0.1471
0.1482
0.1514

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

a)

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

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. 6. Piezoelectric vibration devices with optimum design parameters

Piezoelectric vibration devices  with optimum design parameters

Fig. 7. Measurement 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. 8. Measurement of tip displacement with high-speed imaging instrument

 Measurement of tip displacement with high-speed imaging instrument

Fig. 9. 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

a)

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

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. 10. Measurement 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. 11. Vibration displacement of cutting tool tip

Vibration displacement of cutting tool tip

Fig. 12. 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

a)

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

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 4. Lifetime of piezoelectric bone with and without design optimization

Lifetime (h)
1 < T 30
30 < T 50
50 < T 100
T > 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 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.

References

  1. 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. [Search CrossRef]
  2. 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. [Search CrossRef]
  3. Ying Chen, Rui Kang Design for reliability of ultrasonic bone knife. Advanced Science Letters, Vol. 4, Issue 6-7, 2011, p. 2147-2151. [Search CrossRef]
  4. 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. [Search CrossRef]
  5. Sergio S. Design of an optimized high-power ultrasonic transducer. Muhlen Ultrasonics Symposium, 1990, p. 1631-1633. [Search CrossRef]
  6. Eisner E., Seager J. S. A longitudinally resonant stub for vibrations of large amplitude. Ultrasonics, 1965, p. 88-98. [Search CrossRef]
  7. 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. [Search CrossRef]
  8. Tierce P., Decarpigny J. N., Hamonic B., Decarpigny J. N. Power Sonic and Ultrasonic Transducer Design. Springer, Heidelberg, 1988. [Search CrossRef]
  9. 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. [Search CrossRef]
  10. Han S. P. A globally convergent method for nonlinear programming. Journal of Optimization Theory and Applications, Vol. 22, Issue 3, 1977, p. 297-309. [Search CrossRef]
  11. Powell M. J. D. Algorithms for nolinear constraints that use lagrangian function. Mathematical Programming, Vol. 14, 1980, p. 224-248. [Search CrossRef]
  12. Bartholomew M. Nonlinear optimization with engineering applications. Springer Science Business Media, LLC, 2008. [Search CrossRef]
  13. 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. [Search CrossRef]