ULTRAS 4829
No. of Pages 13, Model 5G
19 May 2014 Ultrasonics xxx (2014) xxx–xxx 1
Contents lists available at ScienceDirect
Ultrasonics journal homepage: www.elsevier.com/locate/ultras 4 5
Torsional ultrasonic transducer computational design optimization
3 6
Q1
7
Dpt. of Structural Mechanics, University of Granada, Politécnico de Fuentenueva, 18071 Granada, Spain
8 9
a r t i c l e
1 5 1 2 12 13 14 15 16 17 18 19 20 21 22 23 24
J. Melchor ⇑, G. Rus
i n f o
Article history: Received 14 August 2013 Received in revised form 29 April 2014 Accepted 1 May 2014 Available online xxxx
Q3
Keywords: Non-destructive evaluation Ultrasound Torsional ultrasound Probability of detection Finite element method Soft tissue mechanics
a b s t r a c t A torsional piezoelectric ultrasonic sensor design is proposed in this paper and computationally tested and optimized to measure shear stiffness properties of soft tissue. These are correlated with a number of pathologies like tumors, hepatic lesions and others. The reason is that, whereas compressibility is predominantly governed by the fluid phase of the tissue, the shear stiffness is dependent on the stroma micro-architecture, which is directly affected by those pathologies. However, diagnostic tools to quantify them are currently not well developed. The first contribution is a new typology of design adapted to quasifluids. A second contribution is the procedure for design optimization, for which an analytical estimate of the Robust Probability Of Detection, called RPOD, is presented for use as optimality criteria. The RPOD is formulated probabilistically to maximize the probability of detecting the least possible pathology while minimizing the effect of noise. The resulting optimal transducer has a resonance frequency of 28 kHz. Ó 2014 Elsevier B.V. All rights reserved.
26 27 28 29 30 31 32 33 34 35 36 37 38
39 40
1. Introduction
41
This work focuses on designing and optimizing a transducer capable of transmitting and receiving torsional ultrasound waves, intended for tissue characterization applications in the field of clinical diagnosis. The application of ultrasonic sensors in the clinical field is mainly covered by compressional waves [1]. However, torsional ultrasonic waves are proposed in this paper for the first time for shear stiffness based medical diagnosis. Torsional ultrasound is chosen as the physical magnitude for several reasons. First, it is a mechanical wave, controlled by and therefore most sensitive to the mechanical parameters than any other indirect measurement. Second, the novelty and impact of this application is rooted in the fact that torsional waves propagate at the S-wave speed cs has been reported [2–4] to be significantly sensitive to consistence changes in quasifluids and soft tissues (p-value < 0.001 measured with ultrasound) and to consistency changes caused by tumors [5,6]. Several authors have reported evidence about the potential of shear dynamic viscoelastic properties to characterize tissue condition [7–10]. We propose a design of a transducer to record information on the mechanical shear modulus in tissues. The reason is that shear stiffness directly correlates with the tissue micro-architecture, which is most sensitive to pathologies, whereas compressional waves are predominantly correlated with the fluid phase of the tissue, which remains mostly constant. The main challenge and novelty resides in the capability of separate P waves and S waves based on the frequency contents of each wave. In contrast to
42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Q2 Q1
⇑ Corresponding author. Tel.: +34 95840037. E-mail addresses: [email protected] (J. Melchor), [email protected] (G. Rus).
commercial shear transducers, torsional movement guarantees no P-wave generation at the boundary of the transducer, with is critical in quasifluids, where cs is at least 10 times lower than cp , and P-waves mask S-waves. This is reinforced by using a radially rigid sensor contact surface. The ensuing design comprises an emitting disk and a receiving ring. Previous torsional transducers exist with resonance frequency between [22–38] kHz, which are of the torsional Bolt-clamped Langevin Type (BLT), and consist of disks of alternating polarity piezoelectric elements and two elastic blocks fitted at each end, in order to generate torsional mechanical vibrations of the desired frequencies with the capability of measuring shear modulus [11,12]. The concept of introducing a radiating head mass resembles the Tonpilz design, but has been modified fundamentally to generate a torsional movement instead of a flexural movement [13]. The transducer design that we propose is composed by a combination of elastic and piezoelectric parts. Piezoelectric materials convert electrical signals into mechanical vibrations and viceversa. Piezoelectric ceramic materials, such as lead zirconate titanate (PZT) are widely used in solid-state actuators and sensors for various areas of nondestructive evaluation and bioengineering such as an optical scanners, precision positioning, noise and vibration sensing and cancellation, and linear motors, [14–19]. For this reason, PZT-5 is adopted for our design. The central focus of our contribution is the optimization of the transducer design with the aim of maximizing the sensitivity to shear properties. There are several optimization criteria of ultrasonic sensors, where the conventional characteristics as impedance of the backing material and the thickness of piezoelectric material have been studied [1,20].
http://dx.doi.org/10.1016/j.ultras.2014.05.001 0041-624X/Ó 2014 Elsevier B.V. All rights reserved.
Q1 Please cite this article in press as: J. Melchor, G. Rus, Torsional ultrasonic transducer computational design optimization, Ultrasonics (2014), http:// dx.doi.org/10.1016/j.ultras.2014.05.001
66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95
ULTRAS 4829
No. of Pages 13, Model 5G
19 May 2014 Q1 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145
146
2
J. Melchor, G. Rus / Ultrasonics xxx (2014) xxx–xxx
Numerical approaches can be used to evaluate transducer designs for each of the design parameters required, to maintain transducer performance within a specified range [21,22]. Our transducer is simulated with the Finite Element Method (FEM) using the FEAP software [23], revealing the main mechanisms, which are in summary a circumferential shear movement of the piezoceramic elements, a torsional relative movement between upper and lower disks and rings, an inertial mass by the disks and rings that control the resonant frequency and movement amplitude, and compressional eigenmodes at frequencies at least 100 times higher than the torsional frequency. The model is subsequently validated with an analytical model, where the simplified differential equations of piezoelectric torsional motions are derived and solved in terms of the circumferential displacement and the electric potential, aimed at unveiling the main correlations between the design parameters of the numerical model, and the response. This strategy has previously been used in different designs [11]. We have introduced dermic and connective tissue in the finite element model in view of the first intended application, but it could be extended or modified to other kinds of soft tissues, whose properties are similar, and the Q4 result too. Our simulated values for the S wave velocity ranged [20, 380] [m/s], and for P-waves [1200, 1800] [m/s], while the thicknesses layers in dermic and connective tissue was allowed to range at [0.3 mm, 0.7 mm]. Genetic Algorithms (GA) [24] and FEM have already been used in the literature to optimize waveform and amplitude spectra of ultrasonic transducers [20,25,26], as well as other optimization criteria [17,26,27]. The issue of the probability of detection (POD) has only been addressed independently, under the name of identifiability, in statistics and mathematics, with a wide application in chemistry and physics. In the field of nondestructive characterization, the concept of POD has been discussed under the name of identifiability in relationship with the number of measurements and the number of degrees of freedom to establish a necessity condition to obtain a minimum criteria between the number of measurements and unknown variables [28]. Also, the inversion theory under a probabilistic formulation and the introduction of probability density functions in the model has been examined (see [29,30]) to explain the robustness of the inversion. Other perspectives of probabilistic formulations have been presented providing a general overview of the statistical inversion theory, whose Bayesian and Kalman filtering exposition for nonstationary inverse problems is particularly useful, for example in the micro-damage detection field [31]. Some of these rational approaches have been scarcely applied experimentally for the case of damage characterization [32,33]. In this paper, a new design optimality criteria called Robust Probability of Detection (RPOD) is proposed in this paper as a methodology approach. It is defined with the goal of maximizing the transducer sensitivity to the tissue properties while minimizing the sensitivity to noise ratio of the multilayered shear elastic constants. 2. Theory
153
The piezoelectricity constitutive equations are described in this section in the charge-stress and electric field-strain form, because they are necessaries to implement piezoelectric properties. The ultrasonic transducer is made of steel and piezoelectric ceramic taking into account two layers of dermic and connective tissue where the waves are propagated. This process is simulated by FEAP software based on finite element method and it is detailed below.
154
2.1. Piezoelectric formulation
155
Piezoelectricity is an interaction between mechanical and electric fields, whereby electric charge that accumulates in certain
147 148 149 150 151 152
156
solid materials in response to applied mechanical stress. It is a coupling between the elasticity and the electricity that makes the concept of transduction possible. Linear piezoelectricity is described analytically within constitutive equation, defined how stress (T), strain (S), charge-density displacement (D), and electric field (E) interact by two constitutive equations. The constitutive laws are in Strain-Charge form:
158 159 160 161 162 163
164
T ¼ CE S þ e T E
ð1Þ
D ¼ e S eS E
ð2Þ
166
where CE is the piezoelectric stiffness matrix, e is the piezoelectric coupling coefficient matrix, eT is its transposed and eS is permittivity coefficient matrix. The equilibrium and compatibility equations are stated as:
167 168 169 170
171
r D ¼ 0; E ¼ r/
ð3Þ
1 2
rS T ¼ 0; S ¼ ðru þ ruT Þ
ð4Þ
173
where u ¼ ðu1 ; u2 ; u3 Þ denotes the displacement vector field, and / is the electric potential or voltage. Finally, the following standard sign criteria is used: the electric field and stress values are considered positive for the direction of polarization of the material and for tractions, respectively.
174
2.2. Finite element formulation
179
The governing discretized equation of motion of the system using is written in the form,
180
175 176 177 178
181
182
€ þ KU ¼ FðtÞ MU 2 a3 u1 6 ua 7 7 6 Ua ¼ 6 2a 7 4 u3 5
ð5Þ
ð6Þ
a
184
/
€ are the displacement and acceleration vectors where U and U including the three components of the displacements uai and the voltage degree of freedom /a associated to node a, respectively, M is the mass matrix, K is the stiffness matrix, and FðtÞ is the time history of the applied load. The global stiffness matrix of the element can be expressed as
K¼
Z
B
ðeÞ
DB dX
ð7Þ
V
where D is the elasticity matrix that transforms effective strains to stresses including the electric field coupling at every point of the domain. It depends on the piezoelectric constitutive equations, and can be represented in matrix form as described below.
2
T 11
3
185 186 187 188 189 190
191 ðeÞT
2
C E11 C E12 C E13
0
0
0
e11 e21
0
32
S11
3
193 194 195 196 197
198
7 6T 7 6C 6 0 0 e12 e22 0 7 7 6 S22 7 6 22 7 6 E12 C E22 C E23 0 7 6 76 7 6 7 6 T 33 7 6 C E13 C E23 C E33 0 6 0 0 e13 e23 0 7 6 S33 7 7 6 7 6 7 6 76 7 6 0 0 C E44 0 0 0 0 e34 7 6 S12 7 6 T 12 7 6 0 7 6 76 7 6 7 6 T 13 7 ¼ 6 0 6 0 0 0 C E55 0 0 0 e35 7 7 6 7 6 S13 7 6 7 6 76 7 6 0 0 0 0 C E66 e16 e26 0 7 6 S23 7 6 T 23 7 6 0 7 6 76 7 6 7 6 D1 7 6 e11 e12 e13 6 0 0 e16 e11 e12 0 7 7 6 7 6 E1 7 6 7 6 76 7 6 0 0 e26 e12 e22 0 5 4 E2 5 4 D2 5 4 e21 e22 e23 0 0 0 e34 e35 D3 0 0 0 e33 E3 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl} D
ð8Þ
Q1 Please cite this article in press as: J. Melchor, G. Rus, Torsional ultrasonic transducer computational design optimization, Ultrasonics (2014), http:// dx.doi.org/10.1016/j.ultras.2014.05.001
157
200
ULTRAS 4829
No. of Pages 13, Model 5G
19 May 2014 Q1 201
and,
202
2
ðeÞ
B
204
Na;1
6 6 0 6 6 0 6 6 a 6N 6 ;2 6 a ¼6 6 N;3 6 6 0 6 6 6 0 6 6 0 4 0
0
0
Na;2
0
0
N a;1
Na;1
0
0
N a;1
Na;3
N a;2
0
0
0
0
0
0
0
Table 1 Preliminar dimensions of the sensor design.
3
7 0 7 7 0 7 7 7 0 7 7 7 0 7 7 7 0 7 7 7 Na;1 7 7 a 7 N;2 5 Na;3
ð9Þ
207
where BðeÞ is the strain–displacement matrix of the element ðeÞ, and the superscript T denotes the transpose operator, being Na;i ; i ¼ 1; 2; 3 the shape functions defined for each type of element.
208
2.3. Implementation
209
The numerical tool selected for solving the response of the model (forward problem) is the Finite Element Method (FEM). A 8-node quadratic finite element is chosen because it improves the ratio between precision of shape functions and the size of the system of equations, with 4 degrees of freedom per node described in Eq. (8). It has been implemented to solve the model given by the previous constitutive equations. The piezoelectric element was developed and implemented in the research academic finite element code FEAP [23] because it is an open code allowing us to introduce our piezoelectric element formulation. The strain-charge to stress-charge transformation has been used to introduce the piezoelectric properties:
205 206
210 211 212 213 214 215 216 217 218 219 220
221
CE ¼ S1 E e¼ d
ð10Þ
S1 E
ð11Þ S1 E
T
223
eS ¼ eT d
224
2.4. Parametric geometry
225
Due to the difficulty that mechanical parameter measurement and detection of soft tissue pathologies may lead to, a torsional S-wave transducer design is proposed in this paper, where the transmitter and the receiver are fit inside a circular crown that can be encapsulated in a clinical diagnostic device. To model it, the partial differential equations described above require the geometry, material properties and boundary conditions to completed the problem definition. The geometry of the specimen comprises a 90 degree sector of a circular crown 90 degrees. The cyclic boundary conditions ensures that the behavior is equivalent to a complete model. The boundary conditions are completed with the excitation, which is defined by a 100 V peak-to peak voltage difference applied between the upper and lower sides of the piezoelectric ceramics, with a spike function based on Heaviside function time dependency:
226 227 228 229 230 231 232 233 234 235 236 237 238 239 240
d
ð12Þ
241 243 244 245 246 247 248 249
3
J. Melchor, G. Rus / Ultrasonics xxx (2014) xxx–xxx
f ¼ AHðtÞHðD tÞ
t D
ð13Þ
where A is the amplitude, H is Heaviside function, D is the distance to origin in the x axes, and t is the time. The output of the model (signal V), consist of a difference of voltage at the piezo-ceramic blocks ends in the receiver ring. The PZT-5A laminate is considered as prismatic block of size Lx [mm], Ly [mm] and Lz [mm], as shown in Table 1 and Fig. 1.
Description of design parameters
Dimensions (mm)
Label
Piezoelectric ceramic width Piezoelectric ceramic length Piezoelectric ceramic thickness Disk piezoelectric ceramic eccentricity Disk radius Ring piezoelectric ceramic eccentricity Ring width Disk and ring thickness
1 1 2 2.5 4.25 7.5 2 8
pw pl pt dpe dr rpe rw drt
The geometry construction in FEAP consists in 86 blocks comprising the transmitter disk pair, receiver crown pair, set of piezoelectric ceramic blocks and two layers of dermic and connective tissue and piezoelectric ceramic, respectively (see Fig. 2). Mesh size and time increment discretization are chosen in Section 5 based on a convergence study (as shown in Figs. 5 and 6). The torque sensor simulation was pre-dimensioned based on the dimensions tested of the simplified analytical model described in the next section, and summarized in the table below. The top and bottom central disks and external rings material will be selected later, within the optimization procedure. Crystal symmetry class: Uniaxial Density: 7750 kg/m3 Relative permittivity:
eT e0
1730 0 0 ¼ 0 1730 0 ; 0 0 1730 12
where e0 ¼ 8:854 10 Flexibility matrix:
251 252 253 254 255 256 257 258 259 260 261 262 263
264
ð14Þ 266
F=m
267 268
16:4 5:74 7:22 0 0 0 5:74 16:4 7:22 0 0 0 7:22 7:22 18:8 0 0 0 SE ¼ 1012 ½m2 =N 0 0 47:5 0 0 0 0 0 0 0 47:5 0 0 0 0 0 0 44:3
269
ð15Þ
271
Piezoelectric coupling:
0 0 0 0 584 0 0 0 584 0 0 1012 ½C=N d¼ 0 171 171 374 0 0 0
272
273
ð16Þ 275 276
3. Simplified model of torsion transducer
277
A simplified analytical model of the fundamental oscillatory movement of the torsion ultrasonic transducer is derived. To this end, a number of assumptions are carried out mainly on the relevance of elements of the design and on their movement shape (eigenmode), in order to arrive at a single degree of freedom system. Second, the piezoelectric element is assumed to have a predominant deformation law uniformly distributed and linearly proportional to the electrical excitation. These assumptions are broken down below to yield the analytical equations of motion:
278
3.1. Assumptions
287
Reduction to a single-degree-of-freedom system, where the eigenproblem reduces to a single frequency and a single mode (steady-state movement shape).
Q1 Please cite this article in press as: J. Melchor, G. Rus, Torsional ultrasonic transducer computational design optimization, Ultrasonics (2014), http:// dx.doi.org/10.1016/j.ultras.2014.05.001
250
279 280 281 282 283 284 285 286
289 288 290 291
ULTRAS 4829
No. of Pages 13, Model 5G
19 May 2014 Q1
4
J. Melchor, G. Rus / Ultrasonics xxx (2014) xxx–xxx
Fig. 1. The geometry of the transducer describes a scheme with the outline of the transmitter and reception elements. Projections are shown above.
_________________ Mises Stress 1.85E+01 0.00E+00 1.00E+03 2.00E+03 3.00E+03 4.00E+03 5.00E+03 6.00E+03 7.00E+03 8.00E+03 9.00E+03 1.00E+04 2.99E+05
Time = 1.35E-04
Fig. 2. Different states of torsional transducer at instants t ¼ 9; 18; 117; 135 ls.
292 293 294
295 297 298
299
300 302
303 304 305
306 308
Movement is assumed to be dominated by torsion rotation h in radians. Dynamic equilibrium of torsional moment:
kh þ I€h ¼ 0
ð17Þ
where k is the stiffness in Nm/rad and I is the inertia. Steady-state solution has the form: 0
h ¼ h sinðxtÞ
ð18Þ
as the transient solution is neglected for computing the eigenvalue. Natural frequency (eigenvalue that fulfills the equilibrium (17)):
x¼
rffiffiffi k I
ð19Þ
in rad/s, or f ¼ 2xp in Hz. Stiffness of piezoceramic due to moment:
309 310
311
k¼
M h
ð20Þ
where M is the applied torsional moment. Since each of the n piezoceramic elements is dominated by shear deformation, and they are located at distance d from the center of rotation
M ¼ ndF ¼ ndabrxz
ð21Þ
where F is the resulting force of the shear stress rxz applied over the area a b of the element.
Q1 Please cite this article in press as: J. Melchor, G. Rus, Torsional ultrasonic transducer computational design optimization, Ultrasonics (2014), http:// dx.doi.org/10.1016/j.ultras.2014.05.001
313 314 315 316 317
318 320 321 322
ULTRAS 4829
No. of Pages 13, Model 5G
19 May 2014 Q1
324 325 326 327
328 330 331 332 333 334
335 337
338 339 340 341 342 343
344 346
347 348
349 351
The element described above deforms due to the shear stress the amount xz , which creates the differential displacement u ¼ l=xz between upper and lower sides separated distance l. This links to the rotation of the circular array of piezoceramic elements,
h ¼ u=d ¼ l=dxz
ð22Þ
The effective length of the piezoelectric elements is reduced as eff l ’ 2l to account for the flexibility of the clamping into the eff mass. Hence, h ¼ l =dxz . The piezoceramic behaves linearly elastically,
rxz ¼ G xz
1.5 1
Signal, v [pV]
323
5
J. Melchor, G. Rus / Ultrasonics xxx (2014) xxx–xxx
0.5
0
−0.5
−1
ð23Þ
−1.5
0
20
40
with modified shear stiffness G adding piezoelectric coupling [11]. The inertia against rotation is dominated by the mass blocks of density q, which are either cylindrical or annular (ring-shaped) of radius r. The inertia of the piezoelectric and other elements is comparatively neglectable. Inertia and mass of cylinder:
I¼
mr 2 ; 2
m ¼ pr 2 hq
ð24Þ
where h is the height of the cylinder (in the axial direction). Inertia and mass of ring:
I ¼ mr 2 ;
m ¼ 2prheq
60
80
100
120
140
Time, t [µ s] Fig. 3. example of measurement validating design.
Table 2 Validation analytic design vs. FEM Model
Frequency (kHz)
Relative error to FEM (%)
Analytic disk Analytic ring FEM
27.977 28.180 28.169
0.680 0.039 –
ð25Þ 386
When the transducer contains both a cylinder mass and a ring mass for each of the transmitting and receiving subsystems, their eigenfrequencies should be matched in order to maximize the combined resonance amplification.
ing a dynamic displacement, and (ii) to measure the response (displacements in the time domain). A finite element method model is used in the forward procedure in previous section. We assume the hypothesis that the dynamic behavior of the tissue in its healthy and pathological states is predictable using a well-calibrated model in the sense of faithful reproduction of real measurements. Then, the measured signal is processed to solve the inverse problem, i.e., to determine the changes in the tissue from its original state. A genetic algorithm search tool [24,34] is used to minimize the discrepancy between the experimental readings and the numerically predicted trial response, by means of a cost functional designed to calibrate for coherent uncertainties and noise, and providing maximal robustness and sensitivity. GA is a heuristic search algorithm that mimics natural selection, where a population of candidate solutions are evaluated in terms of their fitness. An iterative process generates, during a number of generations, child populations through the processes of tournament, crossover and mutation and selection, aimed at converging to the fittest solutions.
368
3.2. Validation
4.1. Cost functional
404
369
The resonance frequency of the complete system for the final optimized design parameters is given in Table 5. It is computed from the FEM model by measuring the time between peaks of cycles (see Fig. 3). It is found that they coincide with the frequencies of the sensor calculated from the analytical simplified model with a relative error less than 1% in both cases, shown in Table 2.
The readings from the sensors are time-domain signals denoted by V for the theoretical or synthetic case, and V x for the experimental case. A reading V 0 in the healthy state of the tissue is defined for calibration, and the measurement to analyze is defined as,
405
352 353 354
355
357 358 359
360
362 363 364 365 366
where e is the thickness of the ring (in the radial direction). Subsystem eigenfrequency in the case of cylinder mass. Combining Eqs. (17)–(24),
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 1 u t nabd G f ¼ 2p p leff hr4 q 2
ð26Þ
Subsystem eigenfrequency in the case of ring mass. Combining Eqs. (17)–(25),
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 1 u t nabd G f ¼ 2p 2pleff her 3 q
ð27Þ
387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403
367
370 371 372 373 374
375
4. Inverse procedure
376
The physical principle to mechanically characterize soft tissue is the following. A physical magnitude is propagated along the medium to be analyzed, which distorts the wave until it is measured at an accessible surface (see Fig. 4). The mechanical parameters responsible for the modification of the wave can be inferred from the measured one under certain circumstances by means of the inverse problem theory discussed below. The problem of nondestructive characterization of mechanical tissue properties is solved by a model-based inverse problem (IP) approach that consists of two steps: (i) to excite the system apply-
377 378 379 380 381 382 383 384 385
U¼
V V
ð28Þ
where the RMS values are defined for a discrete function f in time domain f ðt i Þ at N sampling points as,
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 XN1 RMSðf Þ ¼ f ðt i Þ2 i¼0 N
ð29Þ x
A residual c is defined from the misfit or discrepancy U U between the measurements.
c ¼ ðU UÞ
ð30Þ
The cost functional f or fitness function is defined after a residual vector c of size Ni as the quadratic form,
Q1 Please cite this article in press as: J. Melchor, G. Rus, Torsional ultrasonic transducer computational design optimization, Ultrasonics (2014), http:// dx.doi.org/10.1016/j.ultras.2014.05.001
407 408
409
0
RMSðV 0 Þ
x
406
411 412 413
414 416 417 418
419 421 422 423
ULTRAS 4829
No. of Pages 13, Model 5G
19 May 2014 Q1
6
J. Melchor, G. Rus / Ultrasonics xxx (2014) xxx–xxx
Fig. 4. Simplified system for measuring ultrasound wave distortion through tissue and reconstruct mechanical properties.
Fig. 5. Time step convergence.
424 N
426 427 428 429 430 431 432 433
434 436
f ¼
i 1 2 1 1X jcj ¼ c2 2 2 Ni i¼1 i
ð31Þ
It is useful to define an alternative version of the cost functional denoted as f l , with the property of improving the sensitivity while approaching the optimum, just by introducing a logarithm and a small value to ensure its existence. This definition particularly enhances the convergence speed when the minimization is tackled by with genetic algorithms or other random search algorithms (see Gallego and Rus [35]),
f l ¼ logðf þ Þ
ð32Þ
a reference healthy tissue r, denoted by patpr , and the cost functional that collects the effects of those in a scalar function f, as defined above. Nevertheless, we propose a new criterium of POD associated to multivariate pathologies, to consider the case when a sensor should be sensitive to pathologies for a range of different reference or healthy tissues r and a range of tissue layers p, labeled robust POD (RPOD), and defined using a pessimistic criterion as follows,
RPOD ¼ minPODðpat pr Þ ! signalðpat pr Þ2 PODðpat pr Þ ¼ P >1 noiseðpat Þ2
437
4.2. Probability of detection
438
The following POD definition is defined to provide an idea of the probability that a pathology is positively detected, given a monitored tissue specimen, when some noise and system uncertainty are present. The detection and characterization of pathologies is based on the interpretation of the alterations of the measurements due to the presence of the pathology. Other model uncertainties such as density value, dermal layer thickness variability or sensor manufacturing geometrical tolerances, and system noises may also alter these measurements, but they are considered to fall outside the scope of this paper because the effect of those uncertainties is grouped in the inserted measurements noise. We propose to estimate the POD by the probability that the alteration of the measurement caused by the pathology is larger than that caused by the noise. If we label the alteration on the measurement readings caused by the pathology as the SIGNAL component, and the alteration generated by the noise as NOISE, the former definition can be formulated as [34],
439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454
455
POD ¼ P 457 458 459 460 461
jsignalj2 2
jnoisej
! >1
ð33Þ
Furthermore, three variables are be considered in the problem of maximizing the probability of detection (POD), the level of noise, denoted by r, the severity of the pathologies in terms of shear stiffness and damping alteration in different tissue layers p considering
465 466 467 468 469
470
472
e c ; DG ¼ G G e d g; G is shear modulus, where patpr ¼ fDG ¼ G G r r e correspond to reference healthy tissue, and c and d are refer to G connective and dermic tissue respectively. From the definition of the simulated noise, the dependency of the variation of the measurement with increasing noise is also linear. These two considerations about linearity support the proposal that the measurements on a specimen with noise and with pathology can be expressed as Taylor series expansion centered at the case without noise and without pathology, and neglecting higher order terms (hot) than linear, d
d
X dwi dw wi ðpat pr ; rÞ ¼ wi ð0; 0Þ þ pat pr ð0; 0Þ þ r i ð0; 0Þ þ hot dpat pr drffl{zfflfflfflfflfflfflfflffl} p |fflfflfflfflfflfflffl |fflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl} signal
473 474 475 476 477 478 479 480 481 482
483
noise
ð36Þ
485
where i ¼ 1; . . . ; N i are the measuring points for a general case (one in our design). The first term on the right hand side is the measurement at point i without noise nor pathology. The second term is the alteration of that measurement due to the presence of the pathology only, and is labeled SIGNAL, following the reasoning above. The third term is the alteration of the signal originated by the noise only (NOISE). The second term of the Taylor series depends on the sensitivity of the measurements on the pathology, and can be approximated by finite differences,
486
Q1 Please cite this article in press as: J. Melchor, G. Rus, Torsional ultrasonic transducer computational design optimization, Ultrasonics (2014), http:// dx.doi.org/10.1016/j.ultras.2014.05.001
464
ð35Þ
pr
c
463
ð34Þ
patpr
c
462
487 488 489 490 491 492 493 494 495
ULTRAS 4829
No. of Pages 13, Model 5G
19 May 2014 Q1
7
J. Melchor, G. Rus / Ultrasonics xxx (2014) xxx–xxx
Fig. 6. High incremental time step convergence.
1.5
Signal,v [pV]
1
0.5
0
−0.5
−1 Simulated Experimental (10.00% noise) −1.5
0
20
40
60
80
100
120
140
Time,t [µs] Fig. 9. Simulated measurements for model design, noise 10% Fig. 7. Geometry of the blocks and mesh parameters. Table 3 Materials used in the design.
−8
x 10
1.5
Voltage Disp 2 ring Disp 1 disc Disp 2 disc Disp 3 disc Disp 1 tissue Disp 2 tissue Disp 3 tissue
1
0.5
Material
Young modulus (MPa)
Poisson coefficient
Density (kg/m3)
Steel Aluminium Carbon PMMA Ceramic
210 65 150 3 80
0.3 0.3 0.3 0.3 0.3
7800 2700 1500 1200 1800
0
Table 4 Design parameters and their admissible ranges.
−0.5
Design parameters
Range
Reference (mm)
Label
Piezoelectric ceramic width Piezoelectric ceramic length Piezoelectric ceramic thickness Disk piezoelectric ceramic eccentricity Disk radius Ring piezoelectric ceramic eccentricity Ring width Disk and ring ceramic thickness
[0.75, 2] [0.5, 2] [0.5, 4] [1.5, 3.5]
1 1 2 2.5
pw pl pt dpe
[1.75, 5.75] [5.75, 8.5]
4.25 7.5
dr rpe
[1.5, 2.5] [3,13]
2 8
rw drt
−1
−1.5
−2
0
10
20
30
40
50
60
70
80
90
100
Fig. 8. Simulation of evolution of several control points of the model.
496
dwi ~ pr ; 0Þ ¼ wi ~ pr ; 0Þ ðpat ðpat ;pat pr dpat pr ¼
~ pr þ Dpat pr ; 0Þ wi ðpat ~ pr Dpat pr ; 0Þ wi ðpat 2Dpat pr ð37Þ
498
499 500 501 502
~ pr ! 0 is a small pathology used to guarantee that the where pat FEM captures the perturbations produced at small Dpatpr (since the case patpr ¼ 0 with no pathology needs to be computed with a topologically different mesh), in order to compute
wi;patpr ðp0 ; 0Þ wi;patpr ð0; 0Þ. In addition, a central difference scheme, which yields an error of the order OðDpat2pr Þ, becomes available. Since the noise component is linear by definition, a forward difference scheme is adopted, whose OðDrÞ error is sufficient. The third term of the Taylor series can be directly derivated if the definition of POD is assumed,
dwi ¼ ni RMS wFEM ¼ ni RMS i dr
ð38Þ
Q1 Please cite this article in press as: J. Melchor, G. Rus, Torsional ultrasonic transducer computational design optimization, Ultrasonics (2014), http:// dx.doi.org/10.1016/j.ultras.2014.05.001
503 504 505 506 507 508
509 511
ULTRAS 4829
No. of Pages 13, Model 5G
19 May 2014 Q1 512 513
514
8
J. Melchor, G. Rus / Ultrasonics xxx (2014) xxx–xxx
Pm 1
Eqs. (36) and (38) and the relationship jY i j2 ¼ m bined into (33) to obtain,
0 PODpr ¼ P@
517 518 519
can be com-
1
PNi
2 i¼1 ðwi;patpr ð0; 0ÞÞ i > 1A P r2 RMS2 N1i Ni¼1i n2i 1
pat pr 2 N1
PN i 2 ! RMS2 r2 i¼1 ni PODpr ¼ F Sp
520
ð40Þ 522 523
ð39Þ
Using Monte Carlo techniques and error propagation theory the noise in the measurement points can be concluded to follow a normal distribution [34]. Assuming this distribution, the squared sum of the noise ni is known to follow a Chi-square distribution, PN i 2 2 since i¼1 ni ! vNi (e.g. [36]). The parameter of the Chi-square
528
If the noise generator ni is a random variable, the POD is a probability of the stochastic variable pat2pr , described by the cumulative probability density function F,
distribution is the number of degrees of freedom Ni , which in this case is the number of measurement points. In the case that Ni > 10, the Chi-square distribution can be approximated by a pffiffiffiffiffiffiffiffi Gaussian or normal N distribution v2 ðN i Þ NðN i 2=3; 2N i Þ with pffiffiffiffiffiffiffiffi mean Ni 2=3 and standard deviation 2Ni . This approximation in (40) yields,
0
516
2 i¼1 Y i ,
B PN i 2 C B RMS2 r2 i¼1 ni C C B 2 ¼ PBpat pr > XN C i 2C B ðw ð0; 0ÞÞ A @ i ;pat pr i¼1 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Sp
Fig. 10. Example of sensitivity analysis for rpe and dr parameters with total time 200 ls related to turn displacement disk and radial displacement ring.
Q1 Please cite this article in press as: J. Melchor, G. Rus, Torsional ultrasonic transducer computational design optimization, Ultrasonics (2014), http:// dx.doi.org/10.1016/j.ultras.2014.05.001
524 525 526 527 529 530 531 532 533
ULTRAS 4829
No. of Pages 13, Model 5G
19 May 2014 Q1 534
536
537 538 539 540
541
543
544 545
J. Melchor, G. Rus / Ultrasonics xxx (2014) xxx–xxx
" pffiffiffiffiffiffiffiffi# RMS2 r2 ðNi 2=3Þ RMS2 r2 2Ni 2 pat pr ! N ; Sp Sp
ð41Þ
Rx
Since FðxÞ ¼ 1 f ðyÞdy is the cumulative of the normal probability density function f, whose inverse is x ¼ GðFðxÞÞ, the useful pathology severity to noise ratio patpr =r can be expressed from (41) given a POD level as,
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi pat pr 2Ni RMS2 ðNi 2=3Þ ¼ 1 þ G½PODpr Ni 2=3 Sp r
ð42Þ
Note that the analytical expression (41) is only valid for noise with normal distribution at the measurement points.
9
5. Numerical results
546
5.1. Convergence
547
In order to select the space and time discretization size, a convergence study is carried out, aimed at finding an acceptable tradeoff between computational speed and model error. Two discretization parameters are responsible for convergence: (1) time incremental, and (2) element size, which is inversely proportional to the number of elements. The time increment convergence was obtained evaluating in a range within 12 ns and 3200 ns in finite element model and was described by the following figures, where 25 ns was chosen as the optimal time step. P waves speedily oscillating were captured
548
Fig. 11. Example of sensitivity of connective and dermic tissue with total time 200 ls turn radial displacement disk and voltage signal varying secondary wave speed.
Q1 Please cite this article in press as: J. Melchor, G. Rus, Torsional ultrasonic transducer computational design optimization, Ultrasonics (2014), http:// dx.doi.org/10.1016/j.ultras.2014.05.001
549 550 551 552 553 554 555 556 557
ULTRAS 4829
No. of Pages 13, Model 5G
19 May 2014 Q1
10
J. Melchor, G. Rus / Ultrasonics xxx (2014) xxx–xxx
Cost functional, −log (f)
100
2.5
90 80
1
1.5
1
60
0
G (MPa) d
70
50
1
Cost functional, −log (f)
2
40 30 0
0.5
0
5
10
15
20
25
30
35
40
45
10
1
50
10
0
0
−1
0
1
20 Best value Mean
20
30
40
50
60
0
70
80
90
100
Gc (MPa)
Generation
Fig. 12. Cost functional as a function of the parameters with reference Gd and Gc . The cost function is minima for the model (10, 10) noted by +.
1
POD , noise=1.0%
0.9 0.8 0.7 0.6
Gref=300[Pa] d Gref=300[Pa] c Gref=3000[Pa] d Gref=3000[Pa] c Gref=30000[Pa] d Gref=30000[Pa] c RPOD
0.5 0.4 0.3 0.2 3 10
4
10
5
10
6
10
7
10
Model parameters Gc,d [Pa] Fig. 13. Dependency of the POD on the pathology indicator.
Table 5 Final dimensions of transducer based on POD optimization.
558 559 560 561 562 563 564 565 566 567 568 569 570 571
Design parameters
POD optimization (mm)
Label
Piezoelectric width Piezoelectric length Piezoelectric thickness Disk piezoelectric eccentricity Disk radius Ring piezoelectric eccentricity Ring width Disk-ring thickness Disk-ring material
1.9 0.8 2.8 2.7 5.1 5.9 1.6 4.6 Steel
pw pl pt dpe dr rpe rw drt drm
with 12 ns (see Fig. 6 waves slowly oscillating whereas S-waves were only reproduced with 3200 ns. Therefore, P waves have been flittering-out numerically. The study of convergence of the spacial mesh was obtained from twelve mesh refinement parameters q3 ; n; m; n2 ; q5 ; q2 ; q1 ; q4 ; l; p; l1 and l2 . The final combination of parameters was chosen according to an appropriate error to computational time tradeoff, and the final model consists of 1052 nodes, 186 elements and 32 time increments. The following chart describes the geometry of the parametric distribution of the blocks of elements and their mesh parameters (see Fig. 7). The purpose of the numerical results is to obtain conclusions about which experimental design is better in characterizing mechanical tissue properties. To the latter end, three independent
criteria are evaluated numerically: The effects of the excitation and driving frequency combination on (i) the measurements, (ii) the cost function and (iii) the RPOD, are studied for a set of configurations. The scope is to extract some a priori thumb rules that allow to select those with a more accessible minimum in the cost function, and guarantee satisfactory results for a minimization algorithm.
572
5.2. Simulated measurements
579
Consider the specimen described in previous sections, a sample of resulting measurements is shown in Figs. 8 and 9. The simulated signal incorporates a 10% of noise to avoid the inverse crime. The inverse crime consists in using the same model for simulating the measurements than for reconstructing or inverting the unknowns instead of generating the measurements by a different procedure. This is badly regarded in the inverse problem community since a perfect reconstruction is easily achieved, and robustness cannot be assessed, since any small discrepancies between model and reality are not taken into account, which yields an unreal and optimistic result.
580
5.3. Sensitivity analysis
591
The materials of the disk and ring evaluated in the design were steel, aluminum, carbon fiber, ceramic and PMMA with the following mechanical properties (see Table 3 below). These materials were chosen according to their mechanical properties and the quality of the simulated received signal in terms of amplitude. The piezoelectric properties of PZT-5A are provided by the material manufacturer, and the Rayleigh stiffness damping is estimated from our experimental experience. A sensitivity analysis was performed too by assessing the model dimensions according to the ranges chosen with a convergence criteria shown in Table 4. The first part of the sensitivity study for the model parameters concerns the design geometrical dimensions of the sensor within their admissible ranges. It follows that, as the width of the piezo pw piezoelectric width from 0.75 to 2.00 mm increases, the P-wave amplitude for high frequency and frequency increases too. However, the S-wave amplitude at low frequency both outside and inside the ring is indifferent to this range. The same applies to the length of piezoelectric pl piezoelectric ceramic length: the P-wave amplitude is directly correlated with it, for both high and low frequency components of the response.
592
Q1 Please cite this article in press as: J. Melchor, G. Rus, Torsional ultrasonic transducer computational design optimization, Ultrasonics (2014), http:// dx.doi.org/10.1016/j.ultras.2014.05.001
573 574 575 576 577 578
581 582 583 584 585 586 587 588 589 590
593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612
ULTRAS 4829
No. of Pages 13, Model 5G
19 May 2014 Q1
11
J. Melchor, G. Rus / Ultrasonics xxx (2014) xxx–xxx (p/n)|POD=99.9% [MPa]
(p/n)|POD=99.9% [MPa]
2
800 400
1.8
3.5
Piezo thickness [mm]
1200
1.4 1200
Piezo width [mm]
800
1.6
1.2 1
3
80 0
0 80
80
0
40 0
2.5 2 1.5 12
800
00
00
1
12 800
0.8 0.5
1
1.5
0.5
2
800
3
3.5
4
Piezo length [mm]
4.5
(p/n) |POD=99.9% [MPa] 400
11
5
5.5
6
6.5
7
Disc piezo eccentricity [mm] (p/n) |POD=99.9% [MPa] 2.5
0 80
1200
00 800
8
13
13.5
14
14.5
15
15.5
16
16.5
1000
1.5
17
3
4
5
6
7
8
9
10
12
13
Improvement 17199.5% − design: 1.9 0.8 2.8 5.4 10.2 11.8 1.6 4.6 1
20
Optimal Reference
15 10 5 0 −5 −10 −15 −20
11
Disc−ring thickness [mm]
POD , noise=0.00%
Optimality, −log (p/n|POD=99.9%)
Ring piezo eccentricity [mm]
1000
12.5
1100
00
12
400
800
6 400 11.5 12
400
0
800
7
2 110
Ring width [mm]
16 00
80 0 12
9
1000
0 120
Disc radius [mm]
120 0
10
0.8
0.6
0.4
0.2
Best value Mean
0
5
10
15
20
25
30
35
40
45
50
0 4 10
Generation
5
10
6
10
7
10
Model parameter [Pa] Fig. 14. POD optimization.
Fig. 15. Manufactured sensor.
613 614 615 616 617 618 619 620 621
Thereafter, we analyze the sensitivities regarding design ranges for the sensor geometry. dpe denotes the disk piezoelectric ceramic eccentricity from the center of the sensor, which varies from 1:50 to 3.50 mm. The amplitude of the P-wave amplitude is observed to increase with dpe. However the disk radius dr shows an opposite trend, as the P-wave amplitude decreases when dr increases for frequencies above 1 MHz and below 50 kHz. The same happens to the amplitude of waves at low frequency S. The voltage signal increases with the disk radius distance dr.
The next parameter to be analyzed is the ring piezoelectric ceramic eccentricity rpe. Here, the range varies between 5.75 and 8.5 mm, and no change is observed by varying this design parameter. The following parameter was studied for which sensitivity analysis was rw defined as the width of the ring. The ranges studied varied between 1.50 and 2.50 mm and showed no significant differences in the outcome of the movement amplitude either within the interior of the sensor, the ring, or the recorded voltage signal.
Q1 Please cite this article in press as: J. Melchor, G. Rus, Torsional ultrasonic transducer computational design optimization, Ultrasonics (2014), http:// dx.doi.org/10.1016/j.ultras.2014.05.001
622 623 624 625 626 627 628 629 630
ULTRAS 4829
No. of Pages 13, Model 5G
19 May 2014 Q1 631 632 633 634 635 636 637 638
12
J. Melchor, G. Rus / Ultrasonics xxx (2014) xxx–xxx
Notably, only the S wave amplitude at low frequency in radial displacement decreases with this distance. Finally, the disk radius thickness drt generated a decrease in the P wave amplitude at frequencies above 1 MHz and S-wave (circumferential component of the movement) below 50 kHz. The signal voltage amplitude hardly changes for the studied range of values. Some examples of sensitivity analysis have been calculated and are represented below in Fig. 10.
654
5.3.1. Tissue sensitivity analysis Hereinafter, the calculation of the sensitivity of the model parameters that describe the tissue pathology is developed, both for the dermal tissue layer and for the layer of connective tissue, which were chosen as significant within the model. To validate the model, the ranges of values of S-wave velocity in both dermic tissue and connective tissue, as well as the P-wave speed were verified to agree with literature values. Our simulated values ranged for the S wave velocity in [20, 380] [m/s], and for P-waves [1200, 1800] [m/s], while the thicknesses layers in dermic and connective tissue was allowed to range at [0.3 m/m, 0.7 m/m]. (see [37,38] and Fig. 11). By varying the shear modulus of dermal tissue, we find that the amplitude of P waves of frequency above 1 MHz behaves independently for the allowed range of values. Interestingly, the S-wave amplitude at the ring varies, whereas it does not change at the disk.
655
5.4. Effects on cost functional
656
The shape of the cost function provides another subjective way to evaluate the sensitivity of the numerically predicted signals, based on the following criteria:
639 640 641 642 643 644 645 646 647 648 649 650 651 652 653
657 658 659 660 661 662 663 664 665
The existence of local or global minima affects the convergence of the search algorithm. Steep minima are better than those providing soft valleys, due to algorithm convergence performance. Valleys that present shapes close to circular are considered as an indicator of uncoupled mechanical properties of soft tissue parameters.
666
669
Fig. 12 shows a slice of the multidimensional cost function as functions of the parameters Gd and Gc that are equivalent to C E44 ; C E55 and C E66 (see Eq. 8), for configuration of the model.
670
5.5. Effects on POD
671
681
The aforementioned criteria is aimed at evaluating the local behavior of the cost functional, regardless of the noise effects. However, maximizing the POD enables to find the smallest pathology given the presence of experimental noise, independently of the robustness of the convergence of the search. Moreover, the POD serves as a quantifiable and comparable magnitude that can be used as an objective optimality criterion. Fig. 13 shows an example of the POD estimation for one excitation configuration for increasing pathological values, whereas the dependency of the POD on the pathology extent is illustrated for a fixed noise level that amounts to 10%.
682
5.6. POD optimization
683
Once the effects of the POD are identified, the optimization is calculated to demonstrate how the POD serves as an optimality design criterion. The first four graphs show slices of regions around maximum POD over the 8 parameters of the model. The fifth graph shows the convergence of the genetic algorithms search with a population of 20 individuals and 50 generations of the best value.
667 668
672 673 674 675 676 677 678 679 680
684 685 686 687 688
The last graph shows the optimized parameters for the best POD represented in Table 5 and in Fig. 14. Note that the sensitivity of the sensor is increased 172 times with respect to a reference design. Finally, the central frequency of this design results 28 kHz.
689
6. Conclusions
693
An optimal piezoelectric transducer design was carried out by combining a piezoelectric finite element model and a robust analytical estimate of the probability of detection, called RPOD. This allows to estimate the minimum pathology findable given a proposed sensor design, a layered tissue geometry and noise level on measurements. After validation, the RPOD is used as an optimality criterion to feed the used genetic algorithms. An analytical simplified model is formulated and validated with the finite element model, which is aimed at easily predict trends and design parameter dependencies. The design has the ability to clearly separate P and S wave depending on the frequency and the time when the wave arrives at the receiver. Also, note that the simplified analytical model predicts the frequency response of the sensor with less than 1% error, which is considered validated finite element model. S-waves were only reproduced with 3200 ns, flittering-out P waves numerically. Even so, amplitude is not predicted by the simplified model, only central frequency and attenuation effects are not considered in the simplified model. The optimization improves 172 times the design. However, after application of the probability of detection of pathologies in the tissue there are several local minima therefore requires a global search algorithm such as genetic algorithms. A set of sensitivity tests was performed to validate the robustness of the analytical estimates and verify the feasibility, sensitivity and specificity of the designed transducer by the algorithms above. The piezoelectric sensor was eventually manufactured based on the resulting design parameters at Non-Destructive Evaluation Laboratory of the University of Granada by the END lab team (see Fig. 15).
694
Acknowledgements
723
This research was supported by the Ministry of Education DPI2010-17065 and Junta de Andalucía for projects P11-CTS8089 and GGI3000IDIB. The authors would like to thank N. Bochud, L. Peralta, J. Chiachío, M. Chiachío and A. Gómez for their help and encouragement.
724
References
729
[1] S. Van Kervel, J. Thijssen, A calculation scheme for the optimum design of ultrasonic transducers, Ultrasonics 21 (3) (1983) 134–140. [2] B.L. McFarlin et al., Ultrasound insertion loss of the rat cervix, Am. J. Obstet. Gynecol. 193 (6) (2005) S154. [3] B.L. McFarlin, W.D. O’Brien, M.L. Oelze, J.F. Zachary, R.C. White-Traut, Quantitative ultrasound assessment of the rat cervix, J. Ultras. Med. 25 (8) (2006) 1031–1040. [4] M. Lebertre, F. Ossant, L. Vaillant, S. Diridollou, F. Patat, Spatial variation of acoustic parameters in human skin: an in vitro study between 22 and 45 MHz, Ultras. Med. Biol. 28 (5) (2002) 599–615. [5] F. Sebag, J. Vaillant-Lombard, J. Berbis, V. Griset, J. Henry, P. Petit, C. Oliver, Shear wave elastography: a new ultrasound imaging mode for the differential diagnosis of benign and malignant thyroid nodules, J. Clin. Endocrinol. Metab. 95 (12) (2010) 5281–5288. [6] J.M. Chang, W.K. Moon, N. Cho, A. Yi, H.R. Koo, W. Han, D.-Y. Noh, H.-G. Moon, S.J. Kim, Clinical application of shear wave elastography (SWE) in the diagnosis of benign and malignant breast diseases, Breast Cancer Res. Tr. 129 (1) (2011) 89–97. [7] D. Bader, P. Bowker, Mechanical characteristics of skin and underlying tissues in vivo, Biomaterials 4 (4) (1983) 305–308 (ISSN 0142-961). [8] A.S. Ahuja, Tissue as a Voigt body for propagation of ultrasound, Ultrasonic Imaging 1 (2) (1979) 136–143 (ISSN 0161-734). [9] J. Pereira, J. Mansour, B. Davis, Dynamic measurement of the viscoelastic properties of skin, J. Biomech. 24 (2) (1991) 157–162 (ISSN 0021-929).
Q1 Please cite this article in press as: J. Melchor, G. Rus, Torsional ultrasonic transducer computational design optimization, Ultrasonics (2014), http:// dx.doi.org/10.1016/j.ultras.2014.05.001
690 691 692
695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722
725 726 727 728
730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753
ULTRAS 4829
No. of Pages 13, Model 5G
19 May 2014 Q1 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787
Q5
J. Melchor, G. Rus / Ultrasonics xxx (2014) xxx–xxx [10] T. Bigelow, B. McFarlin, W. O’Brien Jr., M. Oelze, In vivo ultrasonic attenuation slope estimates for detecting cervical ripening in rats: preliminary results, J. Acoust. Soc. Am. 123 (3) (2008) 1794–1800 (cited By (since 1996) 0). [11] J. Kim, O. Kwon, Vibration characteristics of piezoelectric torsional transducers, J. Sound Vib. 264 (2) (2003) 453–473. [12] S. Kaneko, S. Nomoto, H. Yamamori, K. Ohya, Load characteristics of a bolted Langevin torsional transducer, Ultrasonics 34 (2) (1996) 239–241. [13] S.C. Butler, A 2.5 kHz magnetostrictive tonpilz sonar transducer design, J. Acoust. Soc. Am. 109 (5) (2001) 2459. [14] J. Li, P. Liu, H. Ding, W. Cao, Modeling characterization and optimization design for PZT transducer used in Near Field Acoustic Levitation, Sens. Actuat. A: Phys. [15] H. Lin, W. Fang, A rib-reinforced micro torsional mirror driven by electrostatic torque generators, Sens. Actuat. A: Phys. 105 (1) (2003) 1–9. [16] A. Glazounov, Q. Zhang, C. Kim, Torsional actuator based on mechanically amplified shear piezoelectric response, Sens. Actuat. A: Phys. 79 (1) (2000) 22– 30. [17] P. Harkness, A. Cardoni, M. Lucas, A comparison of coupling and degenerating mode in longitudinal-torsional step horns, in: IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control. [18] S. Lin, Study on the Langevin piezoelectric ceramic ultrasonic transducer of longitudinal–flexural composite vibrational mode, Ultrasonics 44 (1) (2006) 109–114. [19] M. Aoyagi, T. Suzuki, Y. Tomikawa, Characteristics of a bolt-clamped torsional vibrator using shear-mode piezo-ceramics sandwiched in the axial direction, Ultrasonics 34 (2) (1996) 219–222. [20] J. Thijssen, W. Verhoef, M. Cloostermans, Optimization of ultrasonic transducers, Ultrasonics 23 (1) (1985) 41–46. [21] H. Kwun, W. Jolly, G. Light, E. Wheeler, Effects of variations in design parameters of ultrasonic transducers on performance characteristics, Ultrasonics 26 (2) (1988) 65–72. [22] C. Dutcher, D. Venerus, Compliance effects on the torsional flow of a viscoelastic fluid, J Non-Newtonian Fluid Mech. 150 (2-3) (2008) 154–161. [23] R. Taylor, FEAP-ein Finite Element Analysis Programm, Ing.-Gemeinschaft Klee & Wrigges, 1987.
13
[24] D. Goldberg, Genetic Algorithms in Search, Optimization and Machine Learning, Addison-Wesley Publ, Reading, Massachussets, etc., 1989. [25] F. Wang, X. Zhao, D. Zhang, Y. Wu, Development of novel ultrasonic transducers for microelectronics packaging, J. Mater. Process. Technol. 209 (3) (2009) 1291–1301. [26] A. Schröder, J. Rautenberg, B. Henning, Evaluation of cost functions for FEA based transducer optimization, Phys. Procedia 3 (1) (2010) 1003–1009. [27] E. Heikkola, M. Laitinen, Model-based optimization of ultrasonic transducers, Ultrason. Sonochem. 12 (1) (2005) 53–57. [28] P.-L. Liu, C.-C. Chen, Parametric identification of truss structures by using transient response, J. Sound Vib. 191 (2) (1996) 273–287. [29] A. Tarantola, Inverse problem theory and methods for model parameter estimation, SIAM (2005). [30] A. Tarantola, B. Valette, Inverse problems = quest for information, J. Geophys 50 (3) (1982) 150–170. [31] H. Lam, C.-T. Ng, M. Veidt, Experimental characterization of multiple cracks in a cantilever beam utilizing transient vibration data following a probabilistic approach, J. Sound Vib. 305 (1) (2007) 34–49. [32] G. Rus, S. Lee, S. Chang, S. Wooh, Optimized damage detection of steel plates from noisy impact test, Int. J. Numer. Methods Eng. 68 (7) (2006) 707–727. [33] G. Rus, R. Palma, J. Pérez-Aparicio, Optimal measurement setup for damage detection in piezoelectric plates, Int. J. Eng. Sci. 47 (4) (2009) 554–572. [34] G. Rus, S.Y. Lee, S.Y. Chang, S.C. Wooh, Optimized damage detection of steel plates from noisy impact test, Int. J. Numer. Methods Eng. 68 (2006) 707–727. [35] R. Gallego, G. Rus, Identification of cracks and cavities using the topological sensitivity boundary integral equation, Comput. Mech. 33 (2004) 154–163. [36] L. Rade, B. Westergren, Mathematics Handbook for Science and Engineering, Springer, 1999. [37] E. Barannik, A. Girnyk, V. Tovstiak, A. Marusenko, S. Emelianov, A. Sarvazyan, Doppler ultrasound detection of shear waves remotely induced in tissue phantoms and tissue in vitro, Ultrasonics 40 (1) (2002) 849–852. [38] M. Orescanin, M.A. Qayyum, K.S. Toohey, M.F. Insana, Complex shear modulus of thermally-damaged liver, in: 2009 IEEE International Ultrasonics Symposium (IUS), IEEE, 2009, pp. 127–130.
788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822
Q1 Please cite this article in press as: J. Melchor, G. Rus, Torsional ultrasonic transducer computational design optimization, Ultrasonics (2014), http:// dx.doi.org/10.1016/j.ultras.2014.05.001
No. of Pages 13, Model 5G
19 May 2014 Ultrasonics xxx (2014) xxx–xxx 1
Contents lists available at ScienceDirect
Ultrasonics journal homepage: www.elsevier.com/locate/ultras 4 5
Torsional ultrasonic transducer computational design optimization
3 6
Q1
7
Dpt. of Structural Mechanics, University of Granada, Politécnico de Fuentenueva, 18071 Granada, Spain
8 9
a r t i c l e
1 5 1 2 12 13 14 15 16 17 18 19 20 21 22 23 24
J. Melchor ⇑, G. Rus
i n f o
Article history: Received 14 August 2013 Received in revised form 29 April 2014 Accepted 1 May 2014 Available online xxxx
Q3
Keywords: Non-destructive evaluation Ultrasound Torsional ultrasound Probability of detection Finite element method Soft tissue mechanics
a b s t r a c t A torsional piezoelectric ultrasonic sensor design is proposed in this paper and computationally tested and optimized to measure shear stiffness properties of soft tissue. These are correlated with a number of pathologies like tumors, hepatic lesions and others. The reason is that, whereas compressibility is predominantly governed by the fluid phase of the tissue, the shear stiffness is dependent on the stroma micro-architecture, which is directly affected by those pathologies. However, diagnostic tools to quantify them are currently not well developed. The first contribution is a new typology of design adapted to quasifluids. A second contribution is the procedure for design optimization, for which an analytical estimate of the Robust Probability Of Detection, called RPOD, is presented for use as optimality criteria. The RPOD is formulated probabilistically to maximize the probability of detecting the least possible pathology while minimizing the effect of noise. The resulting optimal transducer has a resonance frequency of 28 kHz. Ó 2014 Elsevier B.V. All rights reserved.
26 27 28 29 30 31 32 33 34 35 36 37 38
39 40
1. Introduction
41
This work focuses on designing and optimizing a transducer capable of transmitting and receiving torsional ultrasound waves, intended for tissue characterization applications in the field of clinical diagnosis. The application of ultrasonic sensors in the clinical field is mainly covered by compressional waves [1]. However, torsional ultrasonic waves are proposed in this paper for the first time for shear stiffness based medical diagnosis. Torsional ultrasound is chosen as the physical magnitude for several reasons. First, it is a mechanical wave, controlled by and therefore most sensitive to the mechanical parameters than any other indirect measurement. Second, the novelty and impact of this application is rooted in the fact that torsional waves propagate at the S-wave speed cs has been reported [2–4] to be significantly sensitive to consistence changes in quasifluids and soft tissues (p-value < 0.001 measured with ultrasound) and to consistency changes caused by tumors [5,6]. Several authors have reported evidence about the potential of shear dynamic viscoelastic properties to characterize tissue condition [7–10]. We propose a design of a transducer to record information on the mechanical shear modulus in tissues. The reason is that shear stiffness directly correlates with the tissue micro-architecture, which is most sensitive to pathologies, whereas compressional waves are predominantly correlated with the fluid phase of the tissue, which remains mostly constant. The main challenge and novelty resides in the capability of separate P waves and S waves based on the frequency contents of each wave. In contrast to
42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65
Q2 Q1
⇑ Corresponding author. Tel.: +34 95840037. E-mail addresses: [email protected] (J. Melchor), [email protected] (G. Rus).
commercial shear transducers, torsional movement guarantees no P-wave generation at the boundary of the transducer, with is critical in quasifluids, where cs is at least 10 times lower than cp , and P-waves mask S-waves. This is reinforced by using a radially rigid sensor contact surface. The ensuing design comprises an emitting disk and a receiving ring. Previous torsional transducers exist with resonance frequency between [22–38] kHz, which are of the torsional Bolt-clamped Langevin Type (BLT), and consist of disks of alternating polarity piezoelectric elements and two elastic blocks fitted at each end, in order to generate torsional mechanical vibrations of the desired frequencies with the capability of measuring shear modulus [11,12]. The concept of introducing a radiating head mass resembles the Tonpilz design, but has been modified fundamentally to generate a torsional movement instead of a flexural movement [13]. The transducer design that we propose is composed by a combination of elastic and piezoelectric parts. Piezoelectric materials convert electrical signals into mechanical vibrations and viceversa. Piezoelectric ceramic materials, such as lead zirconate titanate (PZT) are widely used in solid-state actuators and sensors for various areas of nondestructive evaluation and bioengineering such as an optical scanners, precision positioning, noise and vibration sensing and cancellation, and linear motors, [14–19]. For this reason, PZT-5 is adopted for our design. The central focus of our contribution is the optimization of the transducer design with the aim of maximizing the sensitivity to shear properties. There are several optimization criteria of ultrasonic sensors, where the conventional characteristics as impedance of the backing material and the thickness of piezoelectric material have been studied [1,20].
http://dx.doi.org/10.1016/j.ultras.2014.05.001 0041-624X/Ó 2014 Elsevier B.V. All rights reserved.
Q1 Please cite this article in press as: J. Melchor, G. Rus, Torsional ultrasonic transducer computational design optimization, Ultrasonics (2014), http:// dx.doi.org/10.1016/j.ultras.2014.05.001
66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95
ULTRAS 4829
No. of Pages 13, Model 5G
19 May 2014 Q1 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145
146
2
J. Melchor, G. Rus / Ultrasonics xxx (2014) xxx–xxx
Numerical approaches can be used to evaluate transducer designs for each of the design parameters required, to maintain transducer performance within a specified range [21,22]. Our transducer is simulated with the Finite Element Method (FEM) using the FEAP software [23], revealing the main mechanisms, which are in summary a circumferential shear movement of the piezoceramic elements, a torsional relative movement between upper and lower disks and rings, an inertial mass by the disks and rings that control the resonant frequency and movement amplitude, and compressional eigenmodes at frequencies at least 100 times higher than the torsional frequency. The model is subsequently validated with an analytical model, where the simplified differential equations of piezoelectric torsional motions are derived and solved in terms of the circumferential displacement and the electric potential, aimed at unveiling the main correlations between the design parameters of the numerical model, and the response. This strategy has previously been used in different designs [11]. We have introduced dermic and connective tissue in the finite element model in view of the first intended application, but it could be extended or modified to other kinds of soft tissues, whose properties are similar, and the Q4 result too. Our simulated values for the S wave velocity ranged [20, 380] [m/s], and for P-waves [1200, 1800] [m/s], while the thicknesses layers in dermic and connective tissue was allowed to range at [0.3 mm, 0.7 mm]. Genetic Algorithms (GA) [24] and FEM have already been used in the literature to optimize waveform and amplitude spectra of ultrasonic transducers [20,25,26], as well as other optimization criteria [17,26,27]. The issue of the probability of detection (POD) has only been addressed independently, under the name of identifiability, in statistics and mathematics, with a wide application in chemistry and physics. In the field of nondestructive characterization, the concept of POD has been discussed under the name of identifiability in relationship with the number of measurements and the number of degrees of freedom to establish a necessity condition to obtain a minimum criteria between the number of measurements and unknown variables [28]. Also, the inversion theory under a probabilistic formulation and the introduction of probability density functions in the model has been examined (see [29,30]) to explain the robustness of the inversion. Other perspectives of probabilistic formulations have been presented providing a general overview of the statistical inversion theory, whose Bayesian and Kalman filtering exposition for nonstationary inverse problems is particularly useful, for example in the micro-damage detection field [31]. Some of these rational approaches have been scarcely applied experimentally for the case of damage characterization [32,33]. In this paper, a new design optimality criteria called Robust Probability of Detection (RPOD) is proposed in this paper as a methodology approach. It is defined with the goal of maximizing the transducer sensitivity to the tissue properties while minimizing the sensitivity to noise ratio of the multilayered shear elastic constants. 2. Theory
153
The piezoelectricity constitutive equations are described in this section in the charge-stress and electric field-strain form, because they are necessaries to implement piezoelectric properties. The ultrasonic transducer is made of steel and piezoelectric ceramic taking into account two layers of dermic and connective tissue where the waves are propagated. This process is simulated by FEAP software based on finite element method and it is detailed below.
154
2.1. Piezoelectric formulation
155
Piezoelectricity is an interaction between mechanical and electric fields, whereby electric charge that accumulates in certain
147 148 149 150 151 152
156
solid materials in response to applied mechanical stress. It is a coupling between the elasticity and the electricity that makes the concept of transduction possible. Linear piezoelectricity is described analytically within constitutive equation, defined how stress (T), strain (S), charge-density displacement (D), and electric field (E) interact by two constitutive equations. The constitutive laws are in Strain-Charge form:
158 159 160 161 162 163
164
T ¼ CE S þ e T E
ð1Þ
D ¼ e S eS E
ð2Þ
166
where CE is the piezoelectric stiffness matrix, e is the piezoelectric coupling coefficient matrix, eT is its transposed and eS is permittivity coefficient matrix. The equilibrium and compatibility equations are stated as:
167 168 169 170
171
r D ¼ 0; E ¼ r/
ð3Þ
1 2
rS T ¼ 0; S ¼ ðru þ ruT Þ
ð4Þ
173
where u ¼ ðu1 ; u2 ; u3 Þ denotes the displacement vector field, and / is the electric potential or voltage. Finally, the following standard sign criteria is used: the electric field and stress values are considered positive for the direction of polarization of the material and for tractions, respectively.
174
2.2. Finite element formulation
179
The governing discretized equation of motion of the system using is written in the form,
180
175 176 177 178
181
182
€ þ KU ¼ FðtÞ MU 2 a3 u1 6 ua 7 7 6 Ua ¼ 6 2a 7 4 u3 5
ð5Þ
ð6Þ
a
184
/
€ are the displacement and acceleration vectors where U and U including the three components of the displacements uai and the voltage degree of freedom /a associated to node a, respectively, M is the mass matrix, K is the stiffness matrix, and FðtÞ is the time history of the applied load. The global stiffness matrix of the element can be expressed as
K¼
Z
B
ðeÞ
DB dX
ð7Þ
V
where D is the elasticity matrix that transforms effective strains to stresses including the electric field coupling at every point of the domain. It depends on the piezoelectric constitutive equations, and can be represented in matrix form as described below.
2
T 11
3
185 186 187 188 189 190
191 ðeÞT
2
C E11 C E12 C E13
0
0
0
e11 e21
0
32
S11
3
193 194 195 196 197
198
7 6T 7 6C 6 0 0 e12 e22 0 7 7 6 S22 7 6 22 7 6 E12 C E22 C E23 0 7 6 76 7 6 7 6 T 33 7 6 C E13 C E23 C E33 0 6 0 0 e13 e23 0 7 6 S33 7 7 6 7 6 7 6 76 7 6 0 0 C E44 0 0 0 0 e34 7 6 S12 7 6 T 12 7 6 0 7 6 76 7 6 7 6 T 13 7 ¼ 6 0 6 0 0 0 C E55 0 0 0 e35 7 7 6 7 6 S13 7 6 7 6 76 7 6 0 0 0 0 C E66 e16 e26 0 7 6 S23 7 6 T 23 7 6 0 7 6 76 7 6 7 6 D1 7 6 e11 e12 e13 6 0 0 e16 e11 e12 0 7 7 6 7 6 E1 7 6 7 6 76 7 6 0 0 e26 e12 e22 0 5 4 E2 5 4 D2 5 4 e21 e22 e23 0 0 0 e34 e35 D3 0 0 0 e33 E3 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl} D
ð8Þ
Q1 Please cite this article in press as: J. Melchor, G. Rus, Torsional ultrasonic transducer computational design optimization, Ultrasonics (2014), http:// dx.doi.org/10.1016/j.ultras.2014.05.001
157
200
ULTRAS 4829
No. of Pages 13, Model 5G
19 May 2014 Q1 201
and,
202
2
ðeÞ
B
204
Na;1
6 6 0 6 6 0 6 6 a 6N 6 ;2 6 a ¼6 6 N;3 6 6 0 6 6 6 0 6 6 0 4 0
0
0
Na;2
0
0
N a;1
Na;1
0
0
N a;1
Na;3
N a;2
0
0
0
0
0
0
0
Table 1 Preliminar dimensions of the sensor design.
3
7 0 7 7 0 7 7 7 0 7 7 7 0 7 7 7 0 7 7 7 Na;1 7 7 a 7 N;2 5 Na;3
ð9Þ
207
where BðeÞ is the strain–displacement matrix of the element ðeÞ, and the superscript T denotes the transpose operator, being Na;i ; i ¼ 1; 2; 3 the shape functions defined for each type of element.
208
2.3. Implementation
209
The numerical tool selected for solving the response of the model (forward problem) is the Finite Element Method (FEM). A 8-node quadratic finite element is chosen because it improves the ratio between precision of shape functions and the size of the system of equations, with 4 degrees of freedom per node described in Eq. (8). It has been implemented to solve the model given by the previous constitutive equations. The piezoelectric element was developed and implemented in the research academic finite element code FEAP [23] because it is an open code allowing us to introduce our piezoelectric element formulation. The strain-charge to stress-charge transformation has been used to introduce the piezoelectric properties:
205 206
210 211 212 213 214 215 216 217 218 219 220
221
CE ¼ S1 E e¼ d
ð10Þ
S1 E
ð11Þ S1 E
T
223
eS ¼ eT d
224
2.4. Parametric geometry
225
Due to the difficulty that mechanical parameter measurement and detection of soft tissue pathologies may lead to, a torsional S-wave transducer design is proposed in this paper, where the transmitter and the receiver are fit inside a circular crown that can be encapsulated in a clinical diagnostic device. To model it, the partial differential equations described above require the geometry, material properties and boundary conditions to completed the problem definition. The geometry of the specimen comprises a 90 degree sector of a circular crown 90 degrees. The cyclic boundary conditions ensures that the behavior is equivalent to a complete model. The boundary conditions are completed with the excitation, which is defined by a 100 V peak-to peak voltage difference applied between the upper and lower sides of the piezoelectric ceramics, with a spike function based on Heaviside function time dependency:
226 227 228 229 230 231 232 233 234 235 236 237 238 239 240
d
ð12Þ
241 243 244 245 246 247 248 249
3
J. Melchor, G. Rus / Ultrasonics xxx (2014) xxx–xxx
f ¼ AHðtÞHðD tÞ
t D
ð13Þ
where A is the amplitude, H is Heaviside function, D is the distance to origin in the x axes, and t is the time. The output of the model (signal V), consist of a difference of voltage at the piezo-ceramic blocks ends in the receiver ring. The PZT-5A laminate is considered as prismatic block of size Lx [mm], Ly [mm] and Lz [mm], as shown in Table 1 and Fig. 1.
Description of design parameters
Dimensions (mm)
Label
Piezoelectric ceramic width Piezoelectric ceramic length Piezoelectric ceramic thickness Disk piezoelectric ceramic eccentricity Disk radius Ring piezoelectric ceramic eccentricity Ring width Disk and ring thickness
1 1 2 2.5 4.25 7.5 2 8
pw pl pt dpe dr rpe rw drt
The geometry construction in FEAP consists in 86 blocks comprising the transmitter disk pair, receiver crown pair, set of piezoelectric ceramic blocks and two layers of dermic and connective tissue and piezoelectric ceramic, respectively (see Fig. 2). Mesh size and time increment discretization are chosen in Section 5 based on a convergence study (as shown in Figs. 5 and 6). The torque sensor simulation was pre-dimensioned based on the dimensions tested of the simplified analytical model described in the next section, and summarized in the table below. The top and bottom central disks and external rings material will be selected later, within the optimization procedure. Crystal symmetry class: Uniaxial Density: 7750 kg/m3 Relative permittivity:
eT e0
1730 0 0 ¼ 0 1730 0 ; 0 0 1730 12
where e0 ¼ 8:854 10 Flexibility matrix:
251 252 253 254 255 256 257 258 259 260 261 262 263
264
ð14Þ 266
F=m
267 268
16:4 5:74 7:22 0 0 0 5:74 16:4 7:22 0 0 0 7:22 7:22 18:8 0 0 0 SE ¼ 1012 ½m2 =N 0 0 47:5 0 0 0 0 0 0 0 47:5 0 0 0 0 0 0 44:3
269
ð15Þ
271
Piezoelectric coupling:
0 0 0 0 584 0 0 0 584 0 0 1012 ½C=N d¼ 0 171 171 374 0 0 0
272
273
ð16Þ 275 276
3. Simplified model of torsion transducer
277
A simplified analytical model of the fundamental oscillatory movement of the torsion ultrasonic transducer is derived. To this end, a number of assumptions are carried out mainly on the relevance of elements of the design and on their movement shape (eigenmode), in order to arrive at a single degree of freedom system. Second, the piezoelectric element is assumed to have a predominant deformation law uniformly distributed and linearly proportional to the electrical excitation. These assumptions are broken down below to yield the analytical equations of motion:
278
3.1. Assumptions
287
Reduction to a single-degree-of-freedom system, where the eigenproblem reduces to a single frequency and a single mode (steady-state movement shape).
Q1 Please cite this article in press as: J. Melchor, G. Rus, Torsional ultrasonic transducer computational design optimization, Ultrasonics (2014), http:// dx.doi.org/10.1016/j.ultras.2014.05.001
250
279 280 281 282 283 284 285 286
289 288 290 291
ULTRAS 4829
No. of Pages 13, Model 5G
19 May 2014 Q1
4
J. Melchor, G. Rus / Ultrasonics xxx (2014) xxx–xxx
Fig. 1. The geometry of the transducer describes a scheme with the outline of the transmitter and reception elements. Projections are shown above.
_________________ Mises Stress 1.85E+01 0.00E+00 1.00E+03 2.00E+03 3.00E+03 4.00E+03 5.00E+03 6.00E+03 7.00E+03 8.00E+03 9.00E+03 1.00E+04 2.99E+05
Time = 1.35E-04
Fig. 2. Different states of torsional transducer at instants t ¼ 9; 18; 117; 135 ls.
292 293 294
295 297 298
299
300 302
303 304 305
306 308
Movement is assumed to be dominated by torsion rotation h in radians. Dynamic equilibrium of torsional moment:
kh þ I€h ¼ 0
ð17Þ
where k is the stiffness in Nm/rad and I is the inertia. Steady-state solution has the form: 0
h ¼ h sinðxtÞ
ð18Þ
as the transient solution is neglected for computing the eigenvalue. Natural frequency (eigenvalue that fulfills the equilibrium (17)):
x¼
rffiffiffi k I
ð19Þ
in rad/s, or f ¼ 2xp in Hz. Stiffness of piezoceramic due to moment:
309 310
311
k¼
M h
ð20Þ
where M is the applied torsional moment. Since each of the n piezoceramic elements is dominated by shear deformation, and they are located at distance d from the center of rotation
M ¼ ndF ¼ ndabrxz
ð21Þ
where F is the resulting force of the shear stress rxz applied over the area a b of the element.
Q1 Please cite this article in press as: J. Melchor, G. Rus, Torsional ultrasonic transducer computational design optimization, Ultrasonics (2014), http:// dx.doi.org/10.1016/j.ultras.2014.05.001
313 314 315 316 317
318 320 321 322
ULTRAS 4829
No. of Pages 13, Model 5G
19 May 2014 Q1
324 325 326 327
328 330 331 332 333 334
335 337
338 339 340 341 342 343
344 346
347 348
349 351
The element described above deforms due to the shear stress the amount xz , which creates the differential displacement u ¼ l=xz between upper and lower sides separated distance l. This links to the rotation of the circular array of piezoceramic elements,
h ¼ u=d ¼ l=dxz
ð22Þ
The effective length of the piezoelectric elements is reduced as eff l ’ 2l to account for the flexibility of the clamping into the eff mass. Hence, h ¼ l =dxz . The piezoceramic behaves linearly elastically,
rxz ¼ G xz
1.5 1
Signal, v [pV]
323
5
J. Melchor, G. Rus / Ultrasonics xxx (2014) xxx–xxx
0.5
0
−0.5
−1
ð23Þ
−1.5
0
20
40
with modified shear stiffness G adding piezoelectric coupling [11]. The inertia against rotation is dominated by the mass blocks of density q, which are either cylindrical or annular (ring-shaped) of radius r. The inertia of the piezoelectric and other elements is comparatively neglectable. Inertia and mass of cylinder:
I¼
mr 2 ; 2
m ¼ pr 2 hq
ð24Þ
where h is the height of the cylinder (in the axial direction). Inertia and mass of ring:
I ¼ mr 2 ;
m ¼ 2prheq
60
80
100
120
140
Time, t [µ s] Fig. 3. example of measurement validating design.
Table 2 Validation analytic design vs. FEM Model
Frequency (kHz)
Relative error to FEM (%)
Analytic disk Analytic ring FEM
27.977 28.180 28.169
0.680 0.039 –
ð25Þ 386
When the transducer contains both a cylinder mass and a ring mass for each of the transmitting and receiving subsystems, their eigenfrequencies should be matched in order to maximize the combined resonance amplification.
ing a dynamic displacement, and (ii) to measure the response (displacements in the time domain). A finite element method model is used in the forward procedure in previous section. We assume the hypothesis that the dynamic behavior of the tissue in its healthy and pathological states is predictable using a well-calibrated model in the sense of faithful reproduction of real measurements. Then, the measured signal is processed to solve the inverse problem, i.e., to determine the changes in the tissue from its original state. A genetic algorithm search tool [24,34] is used to minimize the discrepancy between the experimental readings and the numerically predicted trial response, by means of a cost functional designed to calibrate for coherent uncertainties and noise, and providing maximal robustness and sensitivity. GA is a heuristic search algorithm that mimics natural selection, where a population of candidate solutions are evaluated in terms of their fitness. An iterative process generates, during a number of generations, child populations through the processes of tournament, crossover and mutation and selection, aimed at converging to the fittest solutions.
368
3.2. Validation
4.1. Cost functional
404
369
The resonance frequency of the complete system for the final optimized design parameters is given in Table 5. It is computed from the FEM model by measuring the time between peaks of cycles (see Fig. 3). It is found that they coincide with the frequencies of the sensor calculated from the analytical simplified model with a relative error less than 1% in both cases, shown in Table 2.
The readings from the sensors are time-domain signals denoted by V for the theoretical or synthetic case, and V x for the experimental case. A reading V 0 in the healthy state of the tissue is defined for calibration, and the measurement to analyze is defined as,
405
352 353 354
355
357 358 359
360
362 363 364 365 366
where e is the thickness of the ring (in the radial direction). Subsystem eigenfrequency in the case of cylinder mass. Combining Eqs. (17)–(24),
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 1 u t nabd G f ¼ 2p p leff hr4 q 2
ð26Þ
Subsystem eigenfrequency in the case of ring mass. Combining Eqs. (17)–(25),
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u 2 1 u t nabd G f ¼ 2p 2pleff her 3 q
ð27Þ
387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403
367
370 371 372 373 374
375
4. Inverse procedure
376
The physical principle to mechanically characterize soft tissue is the following. A physical magnitude is propagated along the medium to be analyzed, which distorts the wave until it is measured at an accessible surface (see Fig. 4). The mechanical parameters responsible for the modification of the wave can be inferred from the measured one under certain circumstances by means of the inverse problem theory discussed below. The problem of nondestructive characterization of mechanical tissue properties is solved by a model-based inverse problem (IP) approach that consists of two steps: (i) to excite the system apply-
377 378 379 380 381 382 383 384 385
U¼
V V
ð28Þ
where the RMS values are defined for a discrete function f in time domain f ðt i Þ at N sampling points as,
rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 XN1 RMSðf Þ ¼ f ðt i Þ2 i¼0 N
ð29Þ x
A residual c is defined from the misfit or discrepancy U U between the measurements.
c ¼ ðU UÞ
ð30Þ
The cost functional f or fitness function is defined after a residual vector c of size Ni as the quadratic form,
Q1 Please cite this article in press as: J. Melchor, G. Rus, Torsional ultrasonic transducer computational design optimization, Ultrasonics (2014), http:// dx.doi.org/10.1016/j.ultras.2014.05.001
407 408
409
0
RMSðV 0 Þ
x
406
411 412 413
414 416 417 418
419 421 422 423
ULTRAS 4829
No. of Pages 13, Model 5G
19 May 2014 Q1
6
J. Melchor, G. Rus / Ultrasonics xxx (2014) xxx–xxx
Fig. 4. Simplified system for measuring ultrasound wave distortion through tissue and reconstruct mechanical properties.
Fig. 5. Time step convergence.
424 N
426 427 428 429 430 431 432 433
434 436
f ¼
i 1 2 1 1X jcj ¼ c2 2 2 Ni i¼1 i
ð31Þ
It is useful to define an alternative version of the cost functional denoted as f l , with the property of improving the sensitivity while approaching the optimum, just by introducing a logarithm and a small value to ensure its existence. This definition particularly enhances the convergence speed when the minimization is tackled by with genetic algorithms or other random search algorithms (see Gallego and Rus [35]),
f l ¼ logðf þ Þ
ð32Þ
a reference healthy tissue r, denoted by patpr , and the cost functional that collects the effects of those in a scalar function f, as defined above. Nevertheless, we propose a new criterium of POD associated to multivariate pathologies, to consider the case when a sensor should be sensitive to pathologies for a range of different reference or healthy tissues r and a range of tissue layers p, labeled robust POD (RPOD), and defined using a pessimistic criterion as follows,
RPOD ¼ minPODðpat pr Þ ! signalðpat pr Þ2 PODðpat pr Þ ¼ P >1 noiseðpat Þ2
437
4.2. Probability of detection
438
The following POD definition is defined to provide an idea of the probability that a pathology is positively detected, given a monitored tissue specimen, when some noise and system uncertainty are present. The detection and characterization of pathologies is based on the interpretation of the alterations of the measurements due to the presence of the pathology. Other model uncertainties such as density value, dermal layer thickness variability or sensor manufacturing geometrical tolerances, and system noises may also alter these measurements, but they are considered to fall outside the scope of this paper because the effect of those uncertainties is grouped in the inserted measurements noise. We propose to estimate the POD by the probability that the alteration of the measurement caused by the pathology is larger than that caused by the noise. If we label the alteration on the measurement readings caused by the pathology as the SIGNAL component, and the alteration generated by the noise as NOISE, the former definition can be formulated as [34],
439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454
455
POD ¼ P 457 458 459 460 461
jsignalj2 2
jnoisej
! >1
ð33Þ
Furthermore, three variables are be considered in the problem of maximizing the probability of detection (POD), the level of noise, denoted by r, the severity of the pathologies in terms of shear stiffness and damping alteration in different tissue layers p considering
465 466 467 468 469
470
472
e c ; DG ¼ G G e d g; G is shear modulus, where patpr ¼ fDG ¼ G G r r e correspond to reference healthy tissue, and c and d are refer to G connective and dermic tissue respectively. From the definition of the simulated noise, the dependency of the variation of the measurement with increasing noise is also linear. These two considerations about linearity support the proposal that the measurements on a specimen with noise and with pathology can be expressed as Taylor series expansion centered at the case without noise and without pathology, and neglecting higher order terms (hot) than linear, d
d
X dwi dw wi ðpat pr ; rÞ ¼ wi ð0; 0Þ þ pat pr ð0; 0Þ þ r i ð0; 0Þ þ hot dpat pr drffl{zfflfflfflfflfflfflfflffl} p |fflfflfflfflfflfflffl |fflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflffl} signal
473 474 475 476 477 478 479 480 481 482
483
noise
ð36Þ
485
where i ¼ 1; . . . ; N i are the measuring points for a general case (one in our design). The first term on the right hand side is the measurement at point i without noise nor pathology. The second term is the alteration of that measurement due to the presence of the pathology only, and is labeled SIGNAL, following the reasoning above. The third term is the alteration of the signal originated by the noise only (NOISE). The second term of the Taylor series depends on the sensitivity of the measurements on the pathology, and can be approximated by finite differences,
486
Q1 Please cite this article in press as: J. Melchor, G. Rus, Torsional ultrasonic transducer computational design optimization, Ultrasonics (2014), http:// dx.doi.org/10.1016/j.ultras.2014.05.001
464
ð35Þ
pr
c
463
ð34Þ
patpr
c
462
487 488 489 490 491 492 493 494 495
ULTRAS 4829
No. of Pages 13, Model 5G
19 May 2014 Q1
7
J. Melchor, G. Rus / Ultrasonics xxx (2014) xxx–xxx
Fig. 6. High incremental time step convergence.
1.5
Signal,v [pV]
1
0.5
0
−0.5
−1 Simulated Experimental (10.00% noise) −1.5
0
20
40
60
80
100
120
140
Time,t [µs] Fig. 9. Simulated measurements for model design, noise 10% Fig. 7. Geometry of the blocks and mesh parameters. Table 3 Materials used in the design.
−8
x 10
1.5
Voltage Disp 2 ring Disp 1 disc Disp 2 disc Disp 3 disc Disp 1 tissue Disp 2 tissue Disp 3 tissue
1
0.5
Material
Young modulus (MPa)
Poisson coefficient
Density (kg/m3)
Steel Aluminium Carbon PMMA Ceramic
210 65 150 3 80
0.3 0.3 0.3 0.3 0.3
7800 2700 1500 1200 1800
0
Table 4 Design parameters and their admissible ranges.
−0.5
Design parameters
Range
Reference (mm)
Label
Piezoelectric ceramic width Piezoelectric ceramic length Piezoelectric ceramic thickness Disk piezoelectric ceramic eccentricity Disk radius Ring piezoelectric ceramic eccentricity Ring width Disk and ring ceramic thickness
[0.75, 2] [0.5, 2] [0.5, 4] [1.5, 3.5]
1 1 2 2.5
pw pl pt dpe
[1.75, 5.75] [5.75, 8.5]
4.25 7.5
dr rpe
[1.5, 2.5] [3,13]
2 8
rw drt
−1
−1.5
−2
0
10
20
30
40
50
60
70
80
90
100
Fig. 8. Simulation of evolution of several control points of the model.
496
dwi ~ pr ; 0Þ ¼ wi ~ pr ; 0Þ ðpat ðpat ;pat pr dpat pr ¼
~ pr þ Dpat pr ; 0Þ wi ðpat ~ pr Dpat pr ; 0Þ wi ðpat 2Dpat pr ð37Þ
498
499 500 501 502
~ pr ! 0 is a small pathology used to guarantee that the where pat FEM captures the perturbations produced at small Dpatpr (since the case patpr ¼ 0 with no pathology needs to be computed with a topologically different mesh), in order to compute
wi;patpr ðp0 ; 0Þ wi;patpr ð0; 0Þ. In addition, a central difference scheme, which yields an error of the order OðDpat2pr Þ, becomes available. Since the noise component is linear by definition, a forward difference scheme is adopted, whose OðDrÞ error is sufficient. The third term of the Taylor series can be directly derivated if the definition of POD is assumed,
dwi ¼ ni RMS wFEM ¼ ni RMS i dr
ð38Þ
Q1 Please cite this article in press as: J. Melchor, G. Rus, Torsional ultrasonic transducer computational design optimization, Ultrasonics (2014), http:// dx.doi.org/10.1016/j.ultras.2014.05.001
503 504 505 506 507 508
509 511
ULTRAS 4829
No. of Pages 13, Model 5G
19 May 2014 Q1 512 513
514
8
J. Melchor, G. Rus / Ultrasonics xxx (2014) xxx–xxx
Pm 1
Eqs. (36) and (38) and the relationship jY i j2 ¼ m bined into (33) to obtain,
0 PODpr ¼ P@
517 518 519
can be com-
1
PNi
2 i¼1 ðwi;patpr ð0; 0ÞÞ i > 1A P r2 RMS2 N1i Ni¼1i n2i 1
pat pr 2 N1
PN i 2 ! RMS2 r2 i¼1 ni PODpr ¼ F Sp
520
ð40Þ 522 523
ð39Þ
Using Monte Carlo techniques and error propagation theory the noise in the measurement points can be concluded to follow a normal distribution [34]. Assuming this distribution, the squared sum of the noise ni is known to follow a Chi-square distribution, PN i 2 2 since i¼1 ni ! vNi (e.g. [36]). The parameter of the Chi-square
528
If the noise generator ni is a random variable, the POD is a probability of the stochastic variable pat2pr , described by the cumulative probability density function F,
distribution is the number of degrees of freedom Ni , which in this case is the number of measurement points. In the case that Ni > 10, the Chi-square distribution can be approximated by a pffiffiffiffiffiffiffiffi Gaussian or normal N distribution v2 ðN i Þ NðN i 2=3; 2N i Þ with pffiffiffiffiffiffiffiffi mean Ni 2=3 and standard deviation 2Ni . This approximation in (40) yields,
0
516
2 i¼1 Y i ,
B PN i 2 C B RMS2 r2 i¼1 ni C C B 2 ¼ PBpat pr > XN C i 2C B ðw ð0; 0ÞÞ A @ i ;pat pr i¼1 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Sp
Fig. 10. Example of sensitivity analysis for rpe and dr parameters with total time 200 ls related to turn displacement disk and radial displacement ring.
Q1 Please cite this article in press as: J. Melchor, G. Rus, Torsional ultrasonic transducer computational design optimization, Ultrasonics (2014), http:// dx.doi.org/10.1016/j.ultras.2014.05.001
524 525 526 527 529 530 531 532 533
ULTRAS 4829
No. of Pages 13, Model 5G
19 May 2014 Q1 534
536
537 538 539 540
541
543
544 545
J. Melchor, G. Rus / Ultrasonics xxx (2014) xxx–xxx
" pffiffiffiffiffiffiffiffi# RMS2 r2 ðNi 2=3Þ RMS2 r2 2Ni 2 pat pr ! N ; Sp Sp
ð41Þ
Rx
Since FðxÞ ¼ 1 f ðyÞdy is the cumulative of the normal probability density function f, whose inverse is x ¼ GðFðxÞÞ, the useful pathology severity to noise ratio patpr =r can be expressed from (41) given a POD level as,
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffi pat pr 2Ni RMS2 ðNi 2=3Þ ¼ 1 þ G½PODpr Ni 2=3 Sp r
ð42Þ
Note that the analytical expression (41) is only valid for noise with normal distribution at the measurement points.
9
5. Numerical results
546
5.1. Convergence
547
In order to select the space and time discretization size, a convergence study is carried out, aimed at finding an acceptable tradeoff between computational speed and model error. Two discretization parameters are responsible for convergence: (1) time incremental, and (2) element size, which is inversely proportional to the number of elements. The time increment convergence was obtained evaluating in a range within 12 ns and 3200 ns in finite element model and was described by the following figures, where 25 ns was chosen as the optimal time step. P waves speedily oscillating were captured
548
Fig. 11. Example of sensitivity of connective and dermic tissue with total time 200 ls turn radial displacement disk and voltage signal varying secondary wave speed.
Q1 Please cite this article in press as: J. Melchor, G. Rus, Torsional ultrasonic transducer computational design optimization, Ultrasonics (2014), http:// dx.doi.org/10.1016/j.ultras.2014.05.001
549 550 551 552 553 554 555 556 557
ULTRAS 4829
No. of Pages 13, Model 5G
19 May 2014 Q1
10
J. Melchor, G. Rus / Ultrasonics xxx (2014) xxx–xxx
Cost functional, −log (f)
100
2.5
90 80
1
1.5
1
60
0
G (MPa) d
70
50
1
Cost functional, −log (f)
2
40 30 0
0.5
0
5
10
15
20
25
30
35
40
45
10
1
50
10
0
0
−1
0
1
20 Best value Mean
20
30
40
50
60
0
70
80
90
100
Gc (MPa)
Generation
Fig. 12. Cost functional as a function of the parameters with reference Gd and Gc . The cost function is minima for the model (10, 10) noted by +.
1
POD , noise=1.0%
0.9 0.8 0.7 0.6
Gref=300[Pa] d Gref=300[Pa] c Gref=3000[Pa] d Gref=3000[Pa] c Gref=30000[Pa] d Gref=30000[Pa] c RPOD
0.5 0.4 0.3 0.2 3 10
4
10
5
10
6
10
7
10
Model parameters Gc,d [Pa] Fig. 13. Dependency of the POD on the pathology indicator.
Table 5 Final dimensions of transducer based on POD optimization.
558 559 560 561 562 563 564 565 566 567 568 569 570 571
Design parameters
POD optimization (mm)
Label
Piezoelectric width Piezoelectric length Piezoelectric thickness Disk piezoelectric eccentricity Disk radius Ring piezoelectric eccentricity Ring width Disk-ring thickness Disk-ring material
1.9 0.8 2.8 2.7 5.1 5.9 1.6 4.6 Steel
pw pl pt dpe dr rpe rw drt drm
with 12 ns (see Fig. 6 waves slowly oscillating whereas S-waves were only reproduced with 3200 ns. Therefore, P waves have been flittering-out numerically. The study of convergence of the spacial mesh was obtained from twelve mesh refinement parameters q3 ; n; m; n2 ; q5 ; q2 ; q1 ; q4 ; l; p; l1 and l2 . The final combination of parameters was chosen according to an appropriate error to computational time tradeoff, and the final model consists of 1052 nodes, 186 elements and 32 time increments. The following chart describes the geometry of the parametric distribution of the blocks of elements and their mesh parameters (see Fig. 7). The purpose of the numerical results is to obtain conclusions about which experimental design is better in characterizing mechanical tissue properties. To the latter end, three independent
criteria are evaluated numerically: The effects of the excitation and driving frequency combination on (i) the measurements, (ii) the cost function and (iii) the RPOD, are studied for a set of configurations. The scope is to extract some a priori thumb rules that allow to select those with a more accessible minimum in the cost function, and guarantee satisfactory results for a minimization algorithm.
572
5.2. Simulated measurements
579
Consider the specimen described in previous sections, a sample of resulting measurements is shown in Figs. 8 and 9. The simulated signal incorporates a 10% of noise to avoid the inverse crime. The inverse crime consists in using the same model for simulating the measurements than for reconstructing or inverting the unknowns instead of generating the measurements by a different procedure. This is badly regarded in the inverse problem community since a perfect reconstruction is easily achieved, and robustness cannot be assessed, since any small discrepancies between model and reality are not taken into account, which yields an unreal and optimistic result.
580
5.3. Sensitivity analysis
591
The materials of the disk and ring evaluated in the design were steel, aluminum, carbon fiber, ceramic and PMMA with the following mechanical properties (see Table 3 below). These materials were chosen according to their mechanical properties and the quality of the simulated received signal in terms of amplitude. The piezoelectric properties of PZT-5A are provided by the material manufacturer, and the Rayleigh stiffness damping is estimated from our experimental experience. A sensitivity analysis was performed too by assessing the model dimensions according to the ranges chosen with a convergence criteria shown in Table 4. The first part of the sensitivity study for the model parameters concerns the design geometrical dimensions of the sensor within their admissible ranges. It follows that, as the width of the piezo pw piezoelectric width from 0.75 to 2.00 mm increases, the P-wave amplitude for high frequency and frequency increases too. However, the S-wave amplitude at low frequency both outside and inside the ring is indifferent to this range. The same applies to the length of piezoelectric pl piezoelectric ceramic length: the P-wave amplitude is directly correlated with it, for both high and low frequency components of the response.
592
Q1 Please cite this article in press as: J. Melchor, G. Rus, Torsional ultrasonic transducer computational design optimization, Ultrasonics (2014), http:// dx.doi.org/10.1016/j.ultras.2014.05.001
573 574 575 576 577 578
581 582 583 584 585 586 587 588 589 590
593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612
ULTRAS 4829
No. of Pages 13, Model 5G
19 May 2014 Q1
11
J. Melchor, G. Rus / Ultrasonics xxx (2014) xxx–xxx (p/n)|POD=99.9% [MPa]
(p/n)|POD=99.9% [MPa]
2
800 400
1.8
3.5
Piezo thickness [mm]
1200
1.4 1200
Piezo width [mm]
800
1.6
1.2 1
3
80 0
0 80
80
0
40 0
2.5 2 1.5 12
800
00
00
1
12 800
0.8 0.5
1
1.5
0.5
2
800
3
3.5
4
Piezo length [mm]
4.5
(p/n) |POD=99.9% [MPa] 400
11
5
5.5
6
6.5
7
Disc piezo eccentricity [mm] (p/n) |POD=99.9% [MPa] 2.5
0 80
1200
00 800
8
13
13.5
14
14.5
15
15.5
16
16.5
1000
1.5
17
3
4
5
6
7
8
9
10
12
13
Improvement 17199.5% − design: 1.9 0.8 2.8 5.4 10.2 11.8 1.6 4.6 1
20
Optimal Reference
15 10 5 0 −5 −10 −15 −20
11
Disc−ring thickness [mm]
POD , noise=0.00%
Optimality, −log (p/n|POD=99.9%)
Ring piezo eccentricity [mm]
1000
12.5
1100
00
12
400
800
6 400 11.5 12
400
0
800
7
2 110
Ring width [mm]
16 00
80 0 12
9
1000
0 120
Disc radius [mm]
120 0
10
0.8
0.6
0.4
0.2
Best value Mean
0
5
10
15
20
25
30
35
40
45
50
0 4 10
Generation
5
10
6
10
7
10
Model parameter [Pa] Fig. 14. POD optimization.
Fig. 15. Manufactured sensor.
613 614 615 616 617 618 619 620 621
Thereafter, we analyze the sensitivities regarding design ranges for the sensor geometry. dpe denotes the disk piezoelectric ceramic eccentricity from the center of the sensor, which varies from 1:50 to 3.50 mm. The amplitude of the P-wave amplitude is observed to increase with dpe. However the disk radius dr shows an opposite trend, as the P-wave amplitude decreases when dr increases for frequencies above 1 MHz and below 50 kHz. The same happens to the amplitude of waves at low frequency S. The voltage signal increases with the disk radius distance dr.
The next parameter to be analyzed is the ring piezoelectric ceramic eccentricity rpe. Here, the range varies between 5.75 and 8.5 mm, and no change is observed by varying this design parameter. The following parameter was studied for which sensitivity analysis was rw defined as the width of the ring. The ranges studied varied between 1.50 and 2.50 mm and showed no significant differences in the outcome of the movement amplitude either within the interior of the sensor, the ring, or the recorded voltage signal.
Q1 Please cite this article in press as: J. Melchor, G. Rus, Torsional ultrasonic transducer computational design optimization, Ultrasonics (2014), http:// dx.doi.org/10.1016/j.ultras.2014.05.001
622 623 624 625 626 627 628 629 630
ULTRAS 4829
No. of Pages 13, Model 5G
19 May 2014 Q1 631 632 633 634 635 636 637 638
12
J. Melchor, G. Rus / Ultrasonics xxx (2014) xxx–xxx
Notably, only the S wave amplitude at low frequency in radial displacement decreases with this distance. Finally, the disk radius thickness drt generated a decrease in the P wave amplitude at frequencies above 1 MHz and S-wave (circumferential component of the movement) below 50 kHz. The signal voltage amplitude hardly changes for the studied range of values. Some examples of sensitivity analysis have been calculated and are represented below in Fig. 10.
654
5.3.1. Tissue sensitivity analysis Hereinafter, the calculation of the sensitivity of the model parameters that describe the tissue pathology is developed, both for the dermal tissue layer and for the layer of connective tissue, which were chosen as significant within the model. To validate the model, the ranges of values of S-wave velocity in both dermic tissue and connective tissue, as well as the P-wave speed were verified to agree with literature values. Our simulated values ranged for the S wave velocity in [20, 380] [m/s], and for P-waves [1200, 1800] [m/s], while the thicknesses layers in dermic and connective tissue was allowed to range at [0.3 m/m, 0.7 m/m]. (see [37,38] and Fig. 11). By varying the shear modulus of dermal tissue, we find that the amplitude of P waves of frequency above 1 MHz behaves independently for the allowed range of values. Interestingly, the S-wave amplitude at the ring varies, whereas it does not change at the disk.
655
5.4. Effects on cost functional
656
The shape of the cost function provides another subjective way to evaluate the sensitivity of the numerically predicted signals, based on the following criteria:
639 640 641 642 643 644 645 646 647 648 649 650 651 652 653
657 658 659 660 661 662 663 664 665
The existence of local or global minima affects the convergence of the search algorithm. Steep minima are better than those providing soft valleys, due to algorithm convergence performance. Valleys that present shapes close to circular are considered as an indicator of uncoupled mechanical properties of soft tissue parameters.
666
669
Fig. 12 shows a slice of the multidimensional cost function as functions of the parameters Gd and Gc that are equivalent to C E44 ; C E55 and C E66 (see Eq. 8), for configuration of the model.
670
5.5. Effects on POD
671
681
The aforementioned criteria is aimed at evaluating the local behavior of the cost functional, regardless of the noise effects. However, maximizing the POD enables to find the smallest pathology given the presence of experimental noise, independently of the robustness of the convergence of the search. Moreover, the POD serves as a quantifiable and comparable magnitude that can be used as an objective optimality criterion. Fig. 13 shows an example of the POD estimation for one excitation configuration for increasing pathological values, whereas the dependency of the POD on the pathology extent is illustrated for a fixed noise level that amounts to 10%.
682
5.6. POD optimization
683
Once the effects of the POD are identified, the optimization is calculated to demonstrate how the POD serves as an optimality design criterion. The first four graphs show slices of regions around maximum POD over the 8 parameters of the model. The fifth graph shows the convergence of the genetic algorithms search with a population of 20 individuals and 50 generations of the best value.
667 668
672 673 674 675 676 677 678 679 680
684 685 686 687 688
The last graph shows the optimized parameters for the best POD represented in Table 5 and in Fig. 14. Note that the sensitivity of the sensor is increased 172 times with respect to a reference design. Finally, the central frequency of this design results 28 kHz.
689
6. Conclusions
693
An optimal piezoelectric transducer design was carried out by combining a piezoelectric finite element model and a robust analytical estimate of the probability of detection, called RPOD. This allows to estimate the minimum pathology findable given a proposed sensor design, a layered tissue geometry and noise level on measurements. After validation, the RPOD is used as an optimality criterion to feed the used genetic algorithms. An analytical simplified model is formulated and validated with the finite element model, which is aimed at easily predict trends and design parameter dependencies. The design has the ability to clearly separate P and S wave depending on the frequency and the time when the wave arrives at the receiver. Also, note that the simplified analytical model predicts the frequency response of the sensor with less than 1% error, which is considered validated finite element model. S-waves were only reproduced with 3200 ns, flittering-out P waves numerically. Even so, amplitude is not predicted by the simplified model, only central frequency and attenuation effects are not considered in the simplified model. The optimization improves 172 times the design. However, after application of the probability of detection of pathologies in the tissue there are several local minima therefore requires a global search algorithm such as genetic algorithms. A set of sensitivity tests was performed to validate the robustness of the analytical estimates and verify the feasibility, sensitivity and specificity of the designed transducer by the algorithms above. The piezoelectric sensor was eventually manufactured based on the resulting design parameters at Non-Destructive Evaluation Laboratory of the University of Granada by the END lab team (see Fig. 15).
694
Acknowledgements
723
This research was supported by the Ministry of Education DPI2010-17065 and Junta de Andalucía for projects P11-CTS8089 and GGI3000IDIB. The authors would like to thank N. Bochud, L. Peralta, J. Chiachío, M. Chiachío and A. Gómez for their help and encouragement.
724
References
729
[1] S. Van Kervel, J. Thijssen, A calculation scheme for the optimum design of ultrasonic transducers, Ultrasonics 21 (3) (1983) 134–140. [2] B.L. McFarlin et al., Ultrasound insertion loss of the rat cervix, Am. J. Obstet. Gynecol. 193 (6) (2005) S154. [3] B.L. McFarlin, W.D. O’Brien, M.L. Oelze, J.F. Zachary, R.C. White-Traut, Quantitative ultrasound assessment of the rat cervix, J. Ultras. Med. 25 (8) (2006) 1031–1040. [4] M. Lebertre, F. Ossant, L. Vaillant, S. Diridollou, F. Patat, Spatial variation of acoustic parameters in human skin: an in vitro study between 22 and 45 MHz, Ultras. Med. Biol. 28 (5) (2002) 599–615. [5] F. Sebag, J. Vaillant-Lombard, J. Berbis, V. Griset, J. Henry, P. Petit, C. Oliver, Shear wave elastography: a new ultrasound imaging mode for the differential diagnosis of benign and malignant thyroid nodules, J. Clin. Endocrinol. Metab. 95 (12) (2010) 5281–5288. [6] J.M. Chang, W.K. Moon, N. Cho, A. Yi, H.R. Koo, W. Han, D.-Y. Noh, H.-G. Moon, S.J. Kim, Clinical application of shear wave elastography (SWE) in the diagnosis of benign and malignant breast diseases, Breast Cancer Res. Tr. 129 (1) (2011) 89–97. [7] D. Bader, P. Bowker, Mechanical characteristics of skin and underlying tissues in vivo, Biomaterials 4 (4) (1983) 305–308 (ISSN 0142-961). [8] A.S. Ahuja, Tissue as a Voigt body for propagation of ultrasound, Ultrasonic Imaging 1 (2) (1979) 136–143 (ISSN 0161-734). [9] J. Pereira, J. Mansour, B. Davis, Dynamic measurement of the viscoelastic properties of skin, J. Biomech. 24 (2) (1991) 157–162 (ISSN 0021-929).
Q1 Please cite this article in press as: J. Melchor, G. Rus, Torsional ultrasonic transducer computational design optimization, Ultrasonics (2014), http:// dx.doi.org/10.1016/j.ultras.2014.05.001
690 691 692
695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722
725 726 727 728
730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753
ULTRAS 4829
No. of Pages 13, Model 5G
19 May 2014 Q1 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787
Q5
J. Melchor, G. Rus / Ultrasonics xxx (2014) xxx–xxx [10] T. Bigelow, B. McFarlin, W. O’Brien Jr., M. Oelze, In vivo ultrasonic attenuation slope estimates for detecting cervical ripening in rats: preliminary results, J. Acoust. Soc. Am. 123 (3) (2008) 1794–1800 (cited By (since 1996) 0). [11] J. Kim, O. Kwon, Vibration characteristics of piezoelectric torsional transducers, J. Sound Vib. 264 (2) (2003) 453–473. [12] S. Kaneko, S. Nomoto, H. Yamamori, K. Ohya, Load characteristics of a bolted Langevin torsional transducer, Ultrasonics 34 (2) (1996) 239–241. [13] S.C. Butler, A 2.5 kHz magnetostrictive tonpilz sonar transducer design, J. Acoust. Soc. Am. 109 (5) (2001) 2459. [14] J. Li, P. Liu, H. Ding, W. Cao, Modeling characterization and optimization design for PZT transducer used in Near Field Acoustic Levitation, Sens. Actuat. A: Phys. [15] H. Lin, W. Fang, A rib-reinforced micro torsional mirror driven by electrostatic torque generators, Sens. Actuat. A: Phys. 105 (1) (2003) 1–9. [16] A. Glazounov, Q. Zhang, C. Kim, Torsional actuator based on mechanically amplified shear piezoelectric response, Sens. Actuat. A: Phys. 79 (1) (2000) 22– 30. [17] P. Harkness, A. Cardoni, M. Lucas, A comparison of coupling and degenerating mode in longitudinal-torsional step horns, in: IEEE Transactions on Ultrasonics, Ferroelectrics and Frequency Control. [18] S. Lin, Study on the Langevin piezoelectric ceramic ultrasonic transducer of longitudinal–flexural composite vibrational mode, Ultrasonics 44 (1) (2006) 109–114. [19] M. Aoyagi, T. Suzuki, Y. Tomikawa, Characteristics of a bolt-clamped torsional vibrator using shear-mode piezo-ceramics sandwiched in the axial direction, Ultrasonics 34 (2) (1996) 219–222. [20] J. Thijssen, W. Verhoef, M. Cloostermans, Optimization of ultrasonic transducers, Ultrasonics 23 (1) (1985) 41–46. [21] H. Kwun, W. Jolly, G. Light, E. Wheeler, Effects of variations in design parameters of ultrasonic transducers on performance characteristics, Ultrasonics 26 (2) (1988) 65–72. [22] C. Dutcher, D. Venerus, Compliance effects on the torsional flow of a viscoelastic fluid, J Non-Newtonian Fluid Mech. 150 (2-3) (2008) 154–161. [23] R. Taylor, FEAP-ein Finite Element Analysis Programm, Ing.-Gemeinschaft Klee & Wrigges, 1987.
13
[24] D. Goldberg, Genetic Algorithms in Search, Optimization and Machine Learning, Addison-Wesley Publ, Reading, Massachussets, etc., 1989. [25] F. Wang, X. Zhao, D. Zhang, Y. Wu, Development of novel ultrasonic transducers for microelectronics packaging, J. Mater. Process. Technol. 209 (3) (2009) 1291–1301. [26] A. Schröder, J. Rautenberg, B. Henning, Evaluation of cost functions for FEA based transducer optimization, Phys. Procedia 3 (1) (2010) 1003–1009. [27] E. Heikkola, M. Laitinen, Model-based optimization of ultrasonic transducers, Ultrason. Sonochem. 12 (1) (2005) 53–57. [28] P.-L. Liu, C.-C. Chen, Parametric identification of truss structures by using transient response, J. Sound Vib. 191 (2) (1996) 273–287. [29] A. Tarantola, Inverse problem theory and methods for model parameter estimation, SIAM (2005). [30] A. Tarantola, B. Valette, Inverse problems = quest for information, J. Geophys 50 (3) (1982) 150–170. [31] H. Lam, C.-T. Ng, M. Veidt, Experimental characterization of multiple cracks in a cantilever beam utilizing transient vibration data following a probabilistic approach, J. Sound Vib. 305 (1) (2007) 34–49. [32] G. Rus, S. Lee, S. Chang, S. Wooh, Optimized damage detection of steel plates from noisy impact test, Int. J. Numer. Methods Eng. 68 (7) (2006) 707–727. [33] G. Rus, R. Palma, J. Pérez-Aparicio, Optimal measurement setup for damage detection in piezoelectric plates, Int. J. Eng. Sci. 47 (4) (2009) 554–572. [34] G. Rus, S.Y. Lee, S.Y. Chang, S.C. Wooh, Optimized damage detection of steel plates from noisy impact test, Int. J. Numer. Methods Eng. 68 (2006) 707–727. [35] R. Gallego, G. Rus, Identification of cracks and cavities using the topological sensitivity boundary integral equation, Comput. Mech. 33 (2004) 154–163. [36] L. Rade, B. Westergren, Mathematics Handbook for Science and Engineering, Springer, 1999. [37] E. Barannik, A. Girnyk, V. Tovstiak, A. Marusenko, S. Emelianov, A. Sarvazyan, Doppler ultrasound detection of shear waves remotely induced in tissue phantoms and tissue in vitro, Ultrasonics 40 (1) (2002) 849–852. [38] M. Orescanin, M.A. Qayyum, K.S. Toohey, M.F. Insana, Complex shear modulus of thermally-damaged liver, in: 2009 IEEE International Ultrasonics Symposium (IUS), IEEE, 2009, pp. 127–130.
788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822
Q1 Please cite this article in press as: J. Melchor, G. Rus, Torsional ultrasonic transducer computational design optimization, Ultrasonics (2014), http:// dx.doi.org/10.1016/j.ultras.2014.05.001
