Technická univerzita v Liberci Fakulta mechatroniky a mezioborových inženýrských studií

´ univerzita v Liberci Technicka ´ ch inˇ ´ rsky ´ ch studi´ı Fakulta mechatroniky a mezioborovy zeny

´t dizertac ˇn´ı pra ´ ce Autorefera

´ ln´ı analy ´ za diskretizovane ´ho Spektra ´ho rezona ´ toru modelu piezoelektricke

´ lek Petr Ra ˇervenec 2006 Liberec, c

´ univerzita v Liberci Technicka ´ ch inˇ ´ rsky ´ ch studi´ı Fakulta mechatroniky a mezioborovy zeny

ˇ´ırodove ˇdne ´ inˇ ´ rstv´ı Obor: Pr zeny

´ ln´ı analy ´ za diskretizovane ´ho Spektra ´ho rezona ´ toru modelu piezoelektricke

´ lek Petr Ra

ˇ ˇ´ı Maryˇ skolitel: Prof. Ing. Jir ska, CSc. ´ce: Rozsah pra ˇet stran textu: 76 Poc ˇet tabulek: 2 Poc ˇet obra ´ zk˚ Poc u: 35

Anotace Pˇredkl´adan´a dizertaˇcn´ı pr´ace se zab´ yv´a modelov´an´ım rezonanˇcn´ıch charakteristik piezoelektrick´ ych rezon´ator˚ u. Je uveden fyzik´aln´ı popis piezolektrick´ ych materi´al˚ u za pouˇzit´ı line´arn´ıch piezoelektrick´ ych stavov´ ych rovnic a d´ale je definov´ana u ´loha kmit´an´ı piezoelektrick´ ych rezon´ator˚ u. N´asleduje slab´a formulace u ´lohy a jej´ı diskretizace metodou koneˇcn´ ych prvk˚ u, kter´a vede na zobecnˇenou u ´lohu vlastn´ıch ˇc´ısel s ˇr´ıdk´ ymi a strukturovan´ ymi maticemi. Pro jej´ı ˇreˇsen´ı je pouˇzit implicitnˇe restartovan´ y Arnoldiho algoritmus (IRA). Vyˇreˇsen´ım algebraick´e u ´lohy lze nal´ezt rezonanˇcn´ı frekvence piezoelektrick´eho rezon´atoru. Pomoc´ı koeficient˚ u elektromechanick´e vazby jednotliv´ ych m´od˚ u kmit´an´ı jsou posl´eze vybr´any dominantn´ı m´ody. Implementace modelu zahrnuje vytvoˇren´ı geometrie rezon´atoru, sestaven´ı u ´lohy vlastn´ıch ˇc´ısel, jej´ı ˇreˇsen´ı a n´aslednˇe identifikaci jednotliv´ ych m´od˚ u a jejich roztˇr´ıdˇen´ı dle v´ yznamnosti. K tvorbˇe geometrie a s´ıtˇe je pouˇzit volnˇe dostupn´ y software GMSH, pro ˇreˇsen´ı algebraick´e u ´lohy je pouˇzita implementace metody IRA z volnˇe dostupn´e knihovny ARPACK. Ostatn´ı souˇca´sti programov´e implementace jsou prac´ı autora. Program byl testov´an na pod´eln´em a torzn´ım kmit´an´ı piezoelektrick´e tyˇcinky (s velmi dobrou shodou s analytick´ ym ˇreˇsen´ım) a posl´eze na re´aln´e u ´loze tlouˇst’kovˇe stˇriˇzn´ ych kmit˚ u planparaleln´ıho kˇremenn´eho rezon´atoru. Zde byly nalezeny spr´avn´e dominantn´ı m´ody kmit˚ u, pˇriˇcemˇz odchylka od namˇeˇren´ ych rezonanˇcn´ıch frekvenc´ı ˇcin´ı asi 15%. Relativn´ı odstup jednotlv´ ych rezonanˇcn´ıch frekvenc´ı z˚ ustal dobˇre zachov´an. Motivac´ı pr´ace bylo navrˇzen´ı a implementace kompaktn´ıho softwarov´eho modulu, kter´ y by mohl slouˇzit pˇri procesu n´avrhu a v´ yroby piezoelektrick´ ych prvk˚ u s poˇzadovan´ ymi rezonanˇcn´ımi vlastnostmi. Uvedenou odchylku od namˇeˇren´ ych rezonanˇcn´ch frekvenc´ı lze odstranit kalibrac´ı modelu na danou u ´lohu a, pˇri uv´aˇzen´ı kvalitativn´ı spr´avnosti v´ ysledk˚ u modelu, lze tento v dan´e oblasti prakticky vyuˇz´ıvat.

Abstract The aim of the thesis is the modelling of resonant characteristics of piezoelectric resonators. It contains the physical description of piezoelectric materials, using the linear piezoelectric constitutive equations, and then the problem of oscillation of the piezoelectric resonators is defined. It is followed by the weak formulation of the problem and its discretization by the finite element method, which leads to the generalized eigenvalue problem with sparse structured matrices. For their solving, the implicitly restarted Arnoldi algorithm (IRA) is used. The resonant frequencies are subsequently found by solving this eigenvalue problem and the dominant oscillation modes are then selected according their electromechanical coupling coefficients. The computer implementation of the model consists of building of the resonator’s geometry, of compilation of the eigenvalue problem and its solving, and then of the identification of the particular oscillation modes and their sorting according to their significance. For the creation of the geometry and the mesh, the freely available software GMSH is used. For solving the algebraic problem, the implementation of IRA from the ARPACK library is used. The remaining parts of the computer implementation were made by the author. The model was tested on the longitudinal and torsional oscillation of the piezoelectric beam (with very good agreement with the analytical solution) and then used od the practical problem of the thickness shear vibration of the in-plane parallel quartz resonator. In this problem, the right dominant oscillation modes were found. The relative errors from the measurement were about 15%. The relative distances between the particular dominant modes were well preserved. The motivation of the work was the design and implementation of the compact software module, which would be suitable in the design process and the production of the piezoelectric devices with demanded resonant properties. The above mentioned deviation from the measurement can be eliminated by the calibration of the model to the given task and, with respect to the qualitative rightness of the results, the model can practically used.

Obsah ´ Uvod

6

1 Fyzik´ aln´ı popis u ´ lohy Fyzik´aln´ı popis piezoelektrick´ ych materi´al˚ u. . . . . . . . . . . . . . . . . . . . . . . . Kmit´an´ı piezoelektrick´eho kontinua . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8 8 9

2 Slab´ a formulace a diskretizace u ´ lohy 11 Slab´a formulace . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 Diskretizace u ´lohy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3 Oblasti z´ ajmu 13 Voln´e kmit´an´ı . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 V´ ybˇer dominantn´ıch m´od˚ u . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 4 Numerick´ e metody ˇ reˇ sen´ı algebraick´ eu ´ lohy Krylovovsk´e metody . . . . . . . . . . . . . . . . . . . Arnoldiho faktorizace . . . . . . . . . . . . . . . . . . . Implicitnˇe restartovan´ y Arnoldiho algoritmus . . . . . Zobecnˇen´ y probl´em vlastn´ıch ˇc´ısel . . . . . . . . . . . . Pozn´amky - snaha o redukci velikosti algebraick´e u ´lohy

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

15 15 15 16 16 16

5 Poˇ c´ıtaˇ cov´ a implementace modelu

17

6 Testov´ an´ı modelu Kmit´an´ı kˇremenn´e tyˇcinky . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Pod´eln´e kmity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Torzn´ı kmity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

18 18 18 19

7 Aplikace modelu na re´ aln´ eu ´ loze – kmit´ an´ı planparaleln´ıho kˇ remenn´ eho rezon´ atoru 23 V´ ysledky . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 8 Z´ avˇ er

27

Bibliography

28

5

´ Uvod Piezoelektrick´ e materi´ aly a jejich vyuˇ zit´ı Poˇca´tky z´ajmu o piezoelektrick´e materi´aly sahaj´ı aˇz ke konci 19. stolet´ı, kdy bratˇri Curieov´e objevili pˇr´ım´ y piezoelektrick´ y jev – stlaˇcen´ı krystalu vyvolalo elektrick´ y n´aboj na jeho povrchu. Pozdˇeji objevili tak´e zpˇetn´ y piezoelektrick´ y jev – deformaci krystalu zp˚ usobenou elektrick´ ym polem. Postupnˇe bylo objeveno mnoho materi´al˚ u vykazuj´ıc´ıch piezoelektrick´e vlastnosti, jako napˇr. krystalick´ y kˇremen, PZT keramiky, polymern´ı slouˇceniny apod. Ve 20. letech minul´eho stolet´ı se zaˇcaly piezoelektrick´e materi´aly hojnˇe vyuˇz´ıvat v pr˚ umyslov´e v´ yrobˇe. Dnes se pouˇz´ıvaj´ı napˇr. v mˇeˇr´ıc´ı technice, jako senzory veliˇcin, jako gener´atory vln apod. Podstatnou roli v tˇechto aplikac´ıch hraj´ı piezoelektrick´e rezon´atory. y rezon´ator je tenk´a tyˇcinka ˇci destiˇcka vyroben´a Piezoelectrick´ e rezon´ atory Piezoelectrick´ z piezoelektrick´eho materi´alu, na sv´em povrchu opatˇren´a dvˇema ˇci v´ıce elektrodami (napˇr. [27]). V d˚ usledku harmonick´eho buzen´ı elektrick´ ym polem rezon´ator kmit´a. Parametrem, popisuj´ıc´ım chov´an´ı rezon´atoru, jsou jeho rezonanˇcn´ı frekvence. Ne kaˇzd´ y m´od kmit´an´ı lze vybudit stejnˇe snadno, mnoh´e jsou tlumeny d´ıky materi´alov´ ym a tvarov´ ym vlastnostem rezon´atoru. Je d˚ uleˇzit´e nal´ezt takov´e m´ody kmit˚ u, kter´e lze vybudit snadnˇeji neˇz ostatn´ı. Pracovn´ı oblast rezon´atoru se pak pohybuje v okol´ı dominantn´ı frekvence. Piezoelectrick´e rezon´atory se pouˇz´ıvaj´ı napˇr. jako stabiliz´atory frekvenc´ı elektrick´ ych obvod˚ u, frekvenˇcn´ı filtry, senzory nelektrick´ ych veliˇcin apod.

Motivace dizertaˇ cn´ı pr´ ace Rezonanˇcn´ı frekvence piezoelektrick´ ych rezon´ator˚ u se v praxi (ve v´ yrobˇe) obvykle zjiˇst’uj´ı analytick´ ymi ˇci experiment´aln´ımi metodami. Analytick´e metody jsou ovˇsem plnˇe pouˇziteln´e jen pro jednoduch´e tvary rezon´ator˚ u a z´akladn´ı m´ody kmit˚ u. Nev´ yhodou experiment´aln´ıho testov´an´ı jsou jeho n´aklady. Role matematick´ eho modelov´ an´ı Matematick´e modely, v z´avislosti na jejich komplexnosti, mohou pˇreklenout nev´ yhody jin´ ych metod. Matematick´e modelov´an´ı piezoelektrick´ ych materi´al˚ u se zaˇcalo v hojn´e m´ıˇre vyuˇz´ıvat v posledn´ıch 10 aˇz 15 letech, aˇckoli prvn´ı formulace metody koneˇcn´ ych prvk˚ u (MKP) pro tuto oblast byla publikov´ana jiˇz v roce 1970 Allikem a Hughesem [1]. Zm´ın´ıme nˇekolik publikac´ı, popisuj´ıc´ıch r˚ uzn´e oblasti modelov´an´ı piezoelektrick´ ych materi´al˚ u. Roku 1990 byla publikov´ana obecn´a formulace pro 2D a 3D metodu koneˇcn´ ych prvk˚ u [10]. Publikace [24], [8], [17] (v letech 1991-1996) prezentovaly numerick´e v´ ysledky pro aktivn´ı kontrolu vibrac´ı s vyuˇzit´ım piezoelektrick´ ych materi´al˚ u. Pokusy s vyuˇzit´ım koneˇcn´ ych prvk˚ u vyˇsˇs´ıch ˇr´ad˚ u jsou uvedeny napˇr. v [13] nebo [6]. V souˇcasn´e dobˇe se tak´e zaˇc´ınaj´ı vyuˇz´ıvat skoˇrepinov´e

6

prvky vyuˇz´ıvaj´ıc´ı Mindlinovy teorie (napˇr. [25]). Podrobn´a reˇserˇse literatury je uvedena napˇr. v [15]. Uveden´e publikace se vˇetˇsinou bud’ nevˇenuj´ı popisu numerick´eho ˇreˇsen´ı u ´lohy, nebo vyuˇz´ıvaj´ı nˇekter´ y z komerˇcnˇe dostupn´ ych programov´ ych bal´ık˚ u (napˇr. ANSYS ˇci ABAQUS, kter´e maj´ı implementov´any z´akladn´ı piezoelektrick´e moduly). C´ıl dizertaˇ cn´ı pr´ ace V pr´aci je pops´an MKP model piezoelektrick´eho rezon´atoru. Vych´az´ı se z fyzik´aln´ıho popisu piezoelektrick´eho kontinua. Slab´a formulace a n´asledn´a diskretizace probl´emu vede na rozmˇern´ y line´arn´ı syst´em s ˇr´ıdk´ ymi strukturovan´ ymi maticemi, pomoc´ı kter´ ych je definov´an zobecnˇen´ y probl´em vlastn´ıch ˇc´ısel. Jeho ˇreˇsen´ım lze pak nal´ezt rezonanˇcn´ı frekvence. D˚ uraz je kladen na v´ yˇse uveden´ y zobecnˇen´ y probl´em vlastn´ıch ˇc´ısel. V z´avislosti na parametru diskretizace se m˚ uˇze jednat o velmi rozmˇernou u ´lohu, proto je nutn´e zajistit jej´ı efektivn´ı ˇreˇsen´ı. Je tˇreba zm´ınit, ˇze obvykle nen´ı tˇreba zn´at cel´e spektrum dan´eho probl´emu (vˇsechny rezonanˇcn´ı frekvence), ale zaj´ımaj´ı n´as hlavnˇe frekvence dominantn´ı. Pro ˇreˇsen´ı ˇca´steˇcn´eho probl´emu vlastn´ıch ˇc´ısel pak lze s u ´spˇechem pouˇz´ıt nˇekterou z algebraick´ ych metod zaloˇzenou na krylovovsk´ ych podprostorech. C´ılem pr´ace je, na z´akladˇe diskretizovan´eho modelu piezoelektrick´eho rezon´atoru, navrhnout, implementovat a otestovat efektivn´ı numerick´ y algoritmus a vytvoˇrit ucelen´ y programov´ y prostˇredek pro v´ ypoˇcet dominantn´ıch rezonanˇcn´ıch frekvenc´ı a v´ ysledky porovnat s mˇeˇren´ımi. Takto sestaven´ y programov´ y modul je pouˇziteln´ y pro rezon´atory r˚ uzn´ ych tvar˚ u a jeho v´ ysledky mohou b´ yt pouˇzity pˇri n´avrhu a konstrukci rezon´ator˚ u s poˇzadovan´ ymi vlastnostmi.

Struktura dizertaˇ cn´ı pr´ ace Prvn´ı kapitola obsahuje fyzik´an´ı popis piezoelektrick´ ych materi´al˚ u, kombinuj´ıc´ı pohybov´e rovnice pro kmit´an´ı s line´arn´ımi piezoelektrick´ ymi stavov´ ymi rovnicemi. D´ale jsou zde uvedeny vˇsechny nezbytn´e veliˇciny a pojmy. Druh´a kapitola popisuje matematickou formulaci dan´eho probl´emu, d´ale slabou formulaci probl´emu a jej´ı diskretizaci pomoc´ı MKP, vedouc´ı na ˇreˇsen´ı zobecnˇen´eho probl´emu vlastn´ıch ˇc´ısel. Ve tˇret´ı kapitole je obsaˇzen popis r˚ uzn´ ych u ´loh v modelov´an´ı piezoelektrick´ ych materi´al˚ u s d˚ urazem na u ´lohu voln´eho kmit´an´ı a nalezen´ı dominantn´ıch m´od˚ u kmit´an´ı. ˇ Ctvrt´ a kapitola se vˇenuje popisu Krylovovsk´ ych metod a jejich pouˇzit´ı pro ˇreˇsen´ı rozmˇern´ ych probl´emu vlastn´ıch ˇc´ısel, jmenovitˇe popisu implicitnˇe restartovan´eho Arnoldiho algoritmu. Softwarov´e implementaci modelu je vˇenov´ana p´at´a kapitola, obsahuj´ıc´ı struˇcn´ y popis programovac´ıch prac´ı a jednotliv´ ych softwarov´ ych modul˚ u. ˇ a kapitola obsahuje v´ Sest´ ysledky numerick´eho testov´an´ı modelu na analyticky ˇreˇsiteln´ ych u ´loh´ach kmit´an´ı piezoelektrick´e tyˇcinky. V´ ysledky vykazuj´ı dobrou shodu s analytick´ ym ˇreˇsen´ım. ’ Sedm´a kapitola se vˇenuje aplikaci modelu na re´aln´e u ´loze tlouˇst kovˇe stˇriˇzn´ ych kmit˚ u planparaleln´ıho rezon´atoru. Numerick´e v´ ysledky vykazuj´ı kvalitativn´ı shodu s meˇren´ımi (byly identifikov´any spr´avn´e m´ody jako dominantn´ı), absolutn´ı odchylka od meˇren´ı ˇcin´ı okolo 15%. Relativn´ı vzd´alenosti mezi jednotliv´ ymi dominantn´ımi frekvencemi byly dobˇre zachov´any.

7

Fyzik´ aln´ı popis u ´ lohy V t´eto kapitole jsou pops´any z´akladn´ı fyzik´aln´ı vlastnosti piezoelektrick´ ych materi´al˚ u a formulov´an probl´em kmit´an´ı piezoelektrick´eho rezon´atoru. Podrobnˇejˇs´ı popis piezoelektrick´ ych materi´al˚ u lze nal´ezt napˇr. v [27].

Fyzik´ aln´ı popis piezoelektrick´ ych materi´ al˚ u Piezoelektrick´ e stavov´ e rovnice Piezoelektrick´ y krystal pˇredstavuje strukturu, kde na sobˇe vz´ajemnˇe z´avisej´ı deformace a elektrick´e pole. V line´arn´ı teorii je tento vztah pops´an dvojic´ı stavov´ ych rovnic [22]: zobecnˇen´ ym Hookov´ ym z´akonem Tij = cijkl Skl − dkij Ek ,

i, j = 1, 2, 3

a rovnic´ı pˇr´ım´eho piezoelektrick´eho jevu Dk = dkij Sij + εkj Ej ,

k = 1, 2, 3.

Hook˚ uv z´akon vyjadˇruje vztah mezi symetrick´ ymi tenzory napˇet´ı T, deformace S a vektorem intenzity elektrick´eho pole E. Rovnice pˇr´ım´eho piezoelektrick´eho jevu vyjadˇruje vztah mezi vektorem elektrick´eho posunut´ı D, deformac´ı a elektrick´ ym polem. Plat´ı zn´am´e vztahy · ¸ 1 ∂ u˜i ∂ u˜j ∂ ϕ˜ Sij = + , i, j = 1, 2, 3, , Ei = − 2 ∂xj ∂xi ∂xi ˜ = (˜ kde u u1 , u˜2 , u˜3 )T je vektor posunut´ı a ϕ˜ je elektrick´ y potenci´al1 . Koeficienty cijkl , dkij a εij reprezentuj´ı symetrick´e materi´alov´e tenzory, popsan´e v dalˇs´ım odstavci. Materi´ alov´ e vlastnosti Vlastnosti piezoelektrick´eho rezon´atoru u ´zce z´avisej´ı na materi´alu, z nˇehoˇz je rezon´ator vyroben, a jsou charakterizov´any pomoc´ı v´ yˇse uveden´ ych materi´alov´ ych tenzor˚ u. Stavov´e rovnice jsou line´arn´ımi aproximacemi termodynamick´ ych stavov´ ych rovnic a materi´alov´e tenzory v nich hraj´ı roli materi´alov´ ych konstant. Z podm´ınek termodynamick´e stability pak plyne, ˇze tenzory cijkl a εij musej´ı b´ yt symetrick´e a pozitivnˇe definitn´ı [18]. Tenzor elastick´ ych modul˚ u cijkl je symetrick´ y ve vˇsech sv´ ych indexech a m´a obecnˇe jen 21 nez´avisl´ ych sloˇzek. Tenzor piezoelektrick´ ych koeficient˚ u dijk je symetrick´ y ve druh´ ych dvou indexech a m´a obecnˇe 18 nez´avisl´ ych sloˇzek. Tenzor permitivity εij je symetrick´ y a m´a 6 nez´avisl´ ych sloˇzek. Vˇsechny tenzory lze d´ıky jejich symetrii jednoznaˇcnˇe zapsat pomoc´ı zkr´acen´e indexace jako matice [27]. Symetrie a pozitivn´ı definitnost materi´alov´ ych tenzor˚ u maj´ı pˇr´ım´ y vliv na vlastnosti matic vzeˇsl´ ych z pozdˇejˇs´ı diskretizace u ´lohy. 1

Vlnovka znaˇc´ı z´avislost na ˇcase. Pro harmonick´e kmity pak budou oznaˇcen´ı bez vlnovky vyjadˇrovat amplitudy kmit˚ u.

8

Kmit´ an´ı piezoelektrick´ eho kontinua Kmit´ an´ı a rezonance elastick´ ych materi´ al˚ u Kmit´an´ı piezoelektrick´eho syst´emu je zobecnˇen´ım kmit´an´ı ˇcistˇe elastick´ ych materi´al˚ u, proto zde lze uˇz´ıt stejn´e terminologie a princip˚ u. Struˇcnˇe shrˇ nme obecn´e postupy pˇri v´ ypoˇctu kmit´an´ı elastick´eho syst´emu (podrobnˇeji pops´ano napˇr. v [4]). Je-li syst´em (tˇeleso, mnoˇzina hmotn´ ych bod˚ u apod.) vych´ ylen vnˇejˇs´ım z´asahem z rovnov´aˇzn´eho stavu a pot´e je ponech´an bez dalˇs´ıho zat´ıˇzen´ı, zaˇcne kmitat voln´ ymi kmity. Voln´e kmit´an´ı je kombinac´ı vˇsech vlastn´ıch kmit˚ u dan´eho syst´emu s pˇr´ısluˇsn´ ymi vlastn´ımi frekvencemi. Uvaˇzujme syst´em se soustˇredˇen´ ymi parametry s n stupni volnosti. Voln´e netlumen´e kmit´an´ı je pops´ano pohybovou rovnic´ı Mu ¨ + K u = 0, kde uT = [u1 , ..., un ] je vektor posunut´ı, M je hmotnostn´ı matice a K je elastick´a matice. Uvaˇzujeme-li ˇreˇsen´ı pohybov´e rovnice ve tvaru u(t) = veiωt , ˇreˇsen´ı pohybov´e rovnice je ekvivalentn´ı ˇreˇsen´ı zobecnˇen´eho probl´emu vlastn´ıch ˇc´ısel λ ≡ ω2.

(K − λM) v = 0,

ˇ sen´ım u Reˇ ´lohy je mnoˇzina n vlastn´ıch ˇc´ısel a n vlastn´ıch vektor˚ u obsahuj´ıc´ıch rozloˇzen´ı amplitud kmit˚ u. Kmit´a-li syst´em v d˚ usledku (obvykle harmonick´eho) vnˇejˇs´ıho buzen´ı f (t), hovoˇr´ıme o vyucen´em kmit´an´ı. To je pops´ano pohybovou rovnic´ı Mu ¨ + K u = f (t). Je-li frekvence buzen´ı rovna nˇekter´e z vlastn´ıch frekvenc´ı syst´emu, doch´az´ı k rezonanci a pˇr´ısluˇsn´ y m´od kmit´an´ı je velmi silnˇe vybuzen. Kmit´an´ı elastick´ ych materi´al˚ u se ˇreˇs´ı bud’ analytick´ ymi metodami nebo diskretizac´ı u ´lohy (napˇr. na v´ yˇse uveden´ y syst´em se soustˇredˇen´ ymi parametry). Metoda koneˇcn´ ych prvk˚ u (MKP) je dnes nejpouˇz´ıvanˇejˇs´ı diskretizaˇcn´ı metodou. y rezon´ator s hustotou %, Kmit´ an´ı piezoelektrick´ eho kontinua Uvaˇzujme piezoelektrick´ popsan´ y materi´alov´ ymi tenzory. Oznaˇcme objem rezon´atoru jako Ω a jeho povrch jako Γ. Chov´an´ı rezon´atoru je v urˇcit´em ˇcasov´em intervalu [0, T] pops´ano dvˇema diferenci´aln´ımi rovnicemi: Newtonovou pohybovou rovnic´ı %

∂ 2 u˜i ∂Tij = 2 ∂t ∂xj

i = 1, 2, 3,

a kvazistatickou aproximac´ı Maxwellovy rovnice [12] ∇· D=

∂Dj = 0. ∂xj

Po dosazen´ı ze stavov´ ych rovnic dost´av´ame diferenci´aln´ı rovnice pro posunut´ı a elektrick´ y potenci´al µ ¸ ¶ · ∂ 2 u˜i ∂ 1 ∂ u˜k ∂ u˜l ∂ ϕ˜ % 2 = cijkl + + dkij i = 1, 2, 3, ∂t ∂xj 2 ∂xl ∂xk ∂xk µ ¸ ¶ · ∂ 1 ∂ u˜i ∂ u˜j ∂ ϕ˜ 0= dkij + − εkj . ∂xk 2 ∂xj ∂xi ∂xj 9

se sadou poˇca´teˇcn´ıch a okrajov´ ych podm´ınek: u˜i (., 0) u˜i Tij nj ϕ(., ˜ 0) ϕ˜ Dk n k

= = = = = =

ui , 0, fi , ϕ, ϕD , q,

x ∈ Ω, i = 1, 2, 3, i = 1, 2, 3,

x ∈ Γu , x ∈ Γf ,

x ∈ Γϕ , x ∈ Γq ,

kde Γu ∪ Γf = Γ, Γu ∩ Γf = ∅, Γϕ ∪ Γq = Γ, Γϕ ∩ Γq = ∅. Prav´a strana fi znaˇc´ı mechanick´e buzen´ı vnˇejˇs´ı silou, q znaˇc´ı elektrick´e buzen´ı pomoc´ı pˇriloˇzen´eho n´aboje (v pˇr´ıpadˇe voln´ ych kmit˚ u jsou prav´e strany nulov´e). Tˇemito rovnicemi je definov´ana ´ u ´loha kmit´an´ı piezoelektrick´eho kontinua za dan´ ych okrajov´ ych podm´ınek. Uloha bude d´ale numericky ˇreˇsena pomoc´ı diskretizace MKP. Statick´ a deformace piezoelectrick´ eho continua Kromˇe u ´lohy kmit´an´ı lze v oblasti piezoelektrick´ ych materi´al˚ u ˇreˇsit tak´e statickou u ´lohu zat´ıˇzen´ı piezoelektrick´eho prvku konstantn´ı silou ˇci elektrick´ ym polem. Prv´ y pˇr´ıpad nach´az´ı uplatnˇen´ı v oblasti aktu´ator˚ u, druh´ y v oblasti 2 senzor˚ u. Statick´a u ´loha je speci´aln´ım pˇr´ıpadem u ´lohy obecn´e (je zde ∂∂tu˜2i = 0).

10

Slab´ a formulace a diskretizace u ´ lohy V t´eto kapitole je ops´ano odvozen´ı slab´e formulace a n´asledn´a diskretizace pomoc´ı MKP, vedouc´ı na syst´em obyˇcejn´ ych diferenci´aln´ıch rovnic a n´aslednˇe na zobecnˇen´ y probl´em vlastn´ıch ˇc´ısel. Podrobnˇejˇs´ı v´ yklad t´eto oblasti pˇrin´aˇs´ı napˇr. [16], v pr´aci je uveden struˇcn´ y pˇrehled.

Slab´ a formulace Diskretizace u ´lohy je zaloˇzena na tzv. slab´e formulaci. Idea slab´e formulace je pˇrev´est vztahy popsan´e diferenci´aln´ımi rovnicemi na vztahy integr´aln´ı a poˇzadovat jejich splnˇen´ı v urˇcit´em funkcion´aln´ım smyslu, ˇc´ımˇz omez´ıme poˇzadavky na hladkost hledan´eho ˇreˇsen´ı (z´aroveˇ n toto ˇreˇsen´ı m´a st´ale fyzik´aln´ı smysl). (1)

Prostory funkc´ı Ve slab´e formulaci pouˇz´ıv´ame Sobolevovy prostory W2 (Ω) funkc´ı, jejichˇz zobecnˇen´a derivace je ve druh´e mocninˇe integrabiln´ı na oblasti Ω, a d´ale prostor testovac´ıch (1) funkc´ı V (Ω) tvoˇren´ y funkcemi z W2 (Ω). kter´e maj´ı nulovou stopu na hranici Γ. Integr´ aln´ı rovnosti Slabou formulaci pro prostorovou promˇennou odvod´ıme standardn´ım zp˚ usobem [16], zn´asoben´ım diferenci´aln´ıch rovnost´ı testovac´ımi funkcemi w = (w1 , w2 , w3 ) ∈ [V (Ω)]3 , φ ∈ V (Ω), zintegrov´an´ım pˇres oblast Ω, aplikac´ı Greenovy formule a dosazen´ım okrajov´ ych podm´ınek dost´av´ame integr´aln´ı rovnosti µ 2 µ ¶ ¶ D E ∂ u˜i ∂ ϕ˜ % 2 , wi + (cijkl Skl , Rij )Ω + dkij , Rij = fi , w i , ∂t ∂xk Γf Ω Ω µ µ ¶ ¶ D E ∂φ ∂ ϕ˜ ∂φ djik Sik , − εji , = q, φ . ∂xj Ω ∂xi ∂xj Ω Γq h i ∂w i kde Rij = 12 ∂w + ∂xij . Po slab´em ˇreˇsen´ı pak poˇzadujeme splnˇen´ı tˇechto rovnost´ı pro vˇsechny ∂xj volby testovac´ıch funkc´ı.

Diskretizace u ´ lohy MKP konstruuje koneˇcnˇedimenzion´aln´ı aproximaci slab´eho ˇreˇsen´ı. Pouˇz´ıv´ame standardn´ı Lagrangeovsk´e ˇctyˇrstˇeny, aproximuj´ıc´ı oblast Ω, s line´arn´ımi bazick´ ymi funkcemi. Dosazen´ım aproximac´ı posunut´ı a elektrick´eho potenci´alu do integr´aln´ıch rovnost´ı dost´av´ame po u ´pravˇe syst´em obyˇcejn´ ych diferenci´aln´ıch rovnic, kter´ y lze zapsat v blokov´em tvaru: ¨ + KU + PT Φ = F, MU PU − EΦ = Q.

(1) (2)

K je elastick´a matice, M je hmotnostn´ı matice, P je piezoelektrick´a matice a E je elektrick´a matice. Vektory U a Φ obsahuj´ı hodnoty posunut´ı a potenci´alu v uzlech diskretizace. Prav´a strana obsahuje silov´e a elektrick´e buzen´ı kmit˚ u. 11

Po zaveden´ı Dirichletov´ ych okrajov´ ych podm´ınek pro posunut´ı (upevnˇen´ı rezon´atoru) a elektrick´ y potenci´al (um´ıstˇen´ı elektrod, uzemnˇen´ı) dost´av´ame redukovanou soustavu stejn´eho tvaru, kde jsou bloky K, M a E symetrick´e a pozitivnˇe definitn´ı [11] (cel´a matice syst´emu je symetrick´a, obecnˇe indefinitn´ı). Vˇsechny bloky matice jsou ˇr´ıdk´e.

12

Oblasti z´ ajmu V t´eto kapitole jsou pops´any typy u ´loh, se kter´ ymi se lze setkat v oblasti modelov´an´ı piezoelektrick´ ych materi´al˚ u. D˚ uraz je kladen na u ´lohu voln´eho kmit´an´ı, jeˇz bude pozdˇeji ˇreˇsena, ostatn´ı u ´lohy (tlumen´e kmit´an´ı, statick´ y probl´em, aktivn´ı tlumen´ı vibrac´ı) zde slouˇz´ı sp´ıˇse jako dokumentace moˇzn´ ych oblast´ı, kter´e je moˇzno za souˇcasn´eho stavu ˇreˇsit, a nebudu je v autorefer´atu podrobnˇeji uv´adˇet.

Voln´ e kmit´ an´ı Tˇeˇziˇstˇe popisu chov´an´ı piezoelektrick´eho rezon´atoru leˇz´ı v jeho voln´ ych kmitech. Ty urˇcuj´ı, kdy m˚ uˇze syst´em pod vnˇejˇs´ım buzen´ım dos´ahnout rezonance (a na jak´e frekvenci m´a b´ yt vnˇejˇs´ı buzen´ı). V literatuˇre se rozliˇsuj´ı dva typy voln´eho kmit´an´ı piezoelektrick´ ych rezon´ator˚ u. Prvn´ı z nich je kmit´an´ı se zkratovan´ ymi elektrodami, kdy pro tenk´e vrstvy m˚ uˇzeme pˇredpokl´adat Φ = 0 v cel´em objemu rezon´atoru a probl´em se redukuje na standardn´ı u ´lohu zn´amou z kmit´an´ı ˇcistˇe ¨ + KU = 0, se stejn´ elastick´ ych materi´al˚ u, MU ym postupem ˇreˇsen´ı. Druh´ y typ, obecnˇejˇs´ı a l´epe vystihuj´ıc´ı skuteˇcn´ y stav, je voln´e kmit´an´ı s otevˇren´ ymi elektrodami, kdy ˇreˇs´ıme kompletn´ı probl´em vzeˇsl´ y z diskretizace u ´lohy, ¶ µ ¶ µ ¶ µ U 0 K − ω 2 M PT = , Φ 0 P −E kde ω je u ´hlov´a frekvence voln´eho kmit´an´ı. Rozˇs´ıˇrenou metodou pˇri ˇreˇsen´ı t´eto u ´lohy je tzv. statick´a kondenzace, coˇz je dosazen´ı za elektrick´ y potenci´al z druh´e rovnice do prvn´ı, kdy dost´av´ame u ´lohu podobnou standardn´ı u ´loze T −1 2 kmit´an´ı elastick´ ych materi´al˚ u, ovˇsem s pozmˇenˇenou elastickou matic´ı, (K−P E P)U = ω MU. Pˇri tomto postupu ztr´ac´ıme ˇridkou strukturu p˚ uvodn´ı matice, nav´ıc je nutn´e invertovat submatici E. Proto bylo c´ılem pouˇz´ıt metodu, kter´a by za u ´ˇcelem efektivn´ıho v´ ypoˇctu zachovala p˚ uvodn´ı ´ strukturu matic a nevyˇzadovala jej´ı u ´pravy. Ukolem je tedy ˇreˇsit zobecnˇen´ y probl´em vlastn´ıch ˇc´ısel AX = λBX pro

µ A=

K PT P −E

¶

µ , B=

M 0 0 0

¶

µ ,X =

U Φ

¶ , λ = ω2.

V´ ybˇ er dominantn´ıch m´ od˚ u V praxi je snaha vyr´abˇet rezon´atory, u nichˇz vybran´ y m´od kmit´an´ı bude mnohem snadnˇeji vybuditeln´ y neˇz ostatn´ı (tzv. dominantn´ı m´od), resp. u kter´ ych budou dostateˇcn´e odstupy mezi 13

jednotliv´ ymi dominantn´ımi m´ody. M´ıru vybuditelnosti jednotliv´ ych m´od˚ u kmit˚ u lze vyj´adˇrit napˇr´ıklad pomoc´ı jejich koeficientu elektromechanick´e vazby k, kter´ y je definov´an jako k2 =

2 Em , Est Ed

kde

¢ ¢ 1¡ T ¢ 1¡ T 1¡ T U PΦ , Est = U KU , Ed = Φ EΦ 2 2 2 ˇ ım vyˇsˇs´ı je hodnota koeficientu k, t´ım snadnˇeji jsou vz´ajemn´a, elastick´a a dielektrick´a energie. C´ lze dan´ y m´od vybudit. Celkem snadno pak lze spoˇcten´e m´ody kmit˚ u (ˇreˇsen´ım zobecnˇen´eho probl´emu vlastn´ıch ˇc´ısel) roztˇr´ıdit podle jejich v´ yznamnosti. Em =

14

Numerick´ e metody ˇ reˇ sen´ı algebraick´ e u ´ lohy V t´eto kapitole je pops´ana motivace pouˇzit´ı algebraick´ ych metod zaloˇzen´ ych na Krylovovsk´ ych podprostorech a pops´an implicitnˇe restartovan´ y Arnoldiho algoritmus (IRA), kter´ y je vyuˇzit pro ˇreˇsen´ı uveden´eho zobecnˇen´eho probl´emu vlastn´ıch ˇc´ısel.

Krylovovsk´ e metody Jsou dva moˇzn´e pˇr´ıstupy k ˇreˇsen´ı probl´emu vlastn´ıch ˇc´ısel Ax = λx. V prvn´ım pˇr´ıpadˇe m˚ uˇzeme ˇreˇsit kompletn´ı probl´em vlastn´ıch ˇc´ısel (poˇc´ıtat cel´e spektrum), coˇz obn´aˇs´ı transformaci origin´aln´ıho probl´emu na jednoduˇsˇs´ı tvar, z nˇehoˇz lze pak spektrum snadno z´ıskat (napˇr. QR algoritmus, jehoˇz v´ ysledkem je Schur˚ uv rozklad matic). Tento postup je ovˇsem v praxi nepouˇziteln´ y pro velk´e probl´emy, pˇredevˇs´ım kv˚ uli pamˇet’ov´ ym n´arok˚ um. Pokud nejsme nuceni zjiˇst’ovat vˇsechna vlaˇc´ısla matice, m˚ uˇzeme s u ´spˇechem pouˇz´ıt Krylovovsk´e metody (detailn´ı popis tˇechto metod i n´ıˇze uveden´eho IRA lze nal´ezt napˇr. v [20] ˇci [19]). Idea Krylovovsk´ ych metod je v hled´an´ı aproximace vlastn´ıch vektor˚ u jako prvku tzv. 2 k−1 Krylovovsk´eho podprostoru Kk (A, v) = Span{v1 , Av1 , A v1 , ..., A v1 }, tedy jako line´arn´ı kombinaci z vektor˚ u generuj´ıc´ıch dan´ y podprostor. Krylovovsk´ y podprostor obsahuje nejenom informaci o smˇeru vlastn´ıho vektoru pˇr´ısluˇsej´ıc´ımu maxim´aln´ımu vlastn´ımu ˇc´ıslu (jako to dˇel´a mocninn´a metoda), ale vytv´aˇr´ı i informace o vlastn´ıch vektorech v dalˇs´ıch smˇerech. Hledan´ y vlastn´ı vektor nelze v praxi vyj´adˇrit pˇr´ımo jako line´arn´ı kombinaci krylovovsk´e sekvence, protoˇze vektory z t´eto posloupnosti se pro rostouc´ı index k st´avaj´ı line´arnˇe z´avisl´e. Ke konstrukci b´aze Kk , kter´a by byla ortonorm´aln´ı, slouˇz´ı Arnoldiho faktorizace, algoritmus ortogonalizuj´ıc´ı krylovovskou sekvenci.

Arnoldiho faktorizace Arnoldiho algoritmus konstruuje ortonorm´aln´ı b´azi Kk , Vk = (v1 , v2 , ..., vk ) ∈ Rn×(k+1) , dle sch´ematu AVk = Vk Hk + βk vk+1 eT k, kde Hk je matice v horn´ım Hessenbergovˇe tvaru a ek je jednotkov´ y vektor. Spektrum matice Hk lze nal´ezt uˇz snadnˇeji (napˇr. QR algoritmem). M´ame-li vlastn´ı p´ar (λ, s) matice Hk , pak (λ, Vk s) aproximuje vlastn´ı p´ar matice A. Pro zvyˇsuj´ıc´ı se k pak v ide´aln´ım pˇr´ıpadˇe konverguje k vlastn´ımu p´aru matice A.

15

Implicitnˇ e restartovan´ y Arnoldiho algoritmus Obecnˇe nen´ı zn´amo, jak velk´e k je tˇreba k z´ısk´an´ı dobr´e aproximace. Jsme vˇsak z´aroveˇ n omezeni pamˇet’ov´ ymi moˇznostmi (tj. k nem˚ uˇze b´ yt pˇr´ıliˇs velk´e). Tento probl´em lze pˇreklenout tzv. restartem Arnoldiho algoritmu (podrobnˇeji pops´ano ve vlastn´ı pr´aci nebo napˇr. v [19]). M´ame-li n-dimension´aln´ı probl´em vlastn´ıch ˇc´ısel a zaj´ım´a n´as prvn´ıch k vlastn´ıch ˇc´ısel, zvol´ıme nejprve m co nejvˇetˇs´ı tak, aby se n´am veˇsla do pamˇeti Arnoldiho faktorizace velikosti m. Opakovanˇe prov´ad´ıme m-dimension´aln´ı Arnoldiho faktorizaci tak, ˇze nov´ y startovac´ı vektor v1 vznik´a pomoc´ı tzv. polynomi´aln´ıho filtrov´an´ı, kdy redukujeme sloˇzky vektor˚ u n´aleˇzej´ıc´ı smˇer˚ um k + 1, ..., m, a nov´ y startovac´ı vektor udrˇzujeme v k-dimension´aln´ım invariantn´ım podprostoru, kam se postupnˇe stlaˇcuj´ı informace o hledan´ ych vektorech. K tvorbˇe filtraˇcn´ıho polynomu se pouˇzij´ı ta vlastn´ı ˇc´ısla matice Hm z pˇredchoz´ıho kroku, kter´a n´aleˇzej´ı smˇer˚ um, jeˇz chceme odstranit. Implicitn´ı restart. Je pˇr´ıliˇs n´aroˇcn´e postupovat tak, ˇze nejdˇr´ıve aplikujeme filtraˇcn´ı polynom na tvorbu nov´eho startovac´ıho vektoru a pot´e dokonˇc´ıme celou m-rozmˇernou Arnoldiho faktorizaci. Implicitn´ı restart umoˇzn ˇuje vyj´adˇrit pˇr´ımo Arnoldiho rozklad s nov´ ym startovac´ım vektorem pomoc´ı pˇredchoz´ıho rozkladu za pouˇzit´ı implicitnˇe posunut´eho QR algoritmu (podrobnˇeji pops´ano v pr´aci).

Zobecnˇ en´ y probl´ em vlastn´ıch ˇ c´ısel V´ yˇse uveden´ y IRA lze pouˇz´ıt i pro zobecnˇen´ y probl´em vlastn´ıch ˇc´ısel Ax = λBx. Pouˇzijeme-li shift κ (slouˇz´ı k posunut´ı se do ˇc´asti spektra, jeˇz chceme poˇc´ıtat), pomoc´ı tzv. zobecnˇen´e shiftand-invert transformace lze pˇrev´est zobecnˇen´ y probl´em vlastn´ıch ˇc´ısel na obyˇcejn´ y probl´em vlastn´ıch ˇc´ısel, Cx = µx, kde C = (A − κB)−1 B, µ = (λ − κ)−1 . V Arnoldiho algoritmu nyn´ı pˇribude potˇreba ˇreˇsit v kaˇzd´em kroku soustavu rovnic s matic´ı (A − κB) a vznikl´a b´aze je nyn´ı B-ortonorm´aln´ı.

Pozn´ amky - snaha o redukci velikosti algebraick´ eu ´ lohy Nutnost velk´ ych matic Motivac´ı pro pouˇzit´ı IRA bylo pˇredevˇs´ım sn´ıˇzen´ı pamˇet’ov´ ych n´arok˚ u. Origin´aln´ı algebraick´ y probl´em totiˇz mus´ı b´ yt obvykle velk´ ych dimenz´ı, nebot’ jsme nuceni m´ıt co nejjemnˇejˇs´ı diksretizaci, abychom zajistili dobrou aproximaci slab´eho ˇreˇsen´ı. Nav´ıc, ne kaˇzd´ y m´od kmit˚ u lze zachytit libovolnou diskretizac´ı. Matematick´ y model se chov´a jako frekvenˇcn´ı filtr – v z´avislosti na jemnosti diskretizace je schopen postihnout jen m´ody urˇcit´e sloˇzitosti. Pro m´ody komplikovan´ ych tvar˚ u je nutn´a velmi jemn´a s´ıt’, aby tento m´od zachytila (nen´ı zn´am obecn´ y postup, ale v praxi se doporuˇcuje m´ıt alespoˇ n 3 prvky na jednu vlnovou d´elku dan´eho kmitu). Redukce smˇ er˚ u. V pˇr´ıpadˇe, ˇze chceme naj´ıt kmity v dan´em smˇeru ˇci rovinˇe, lze sn´ıˇzit velikost algebraick´e u ´lohy a ˇreˇsit u ´lohu redukovanou na zvolen´ y smˇer ˇci rovinu (pˇri zachov´an´ı informace o vz´ajemn´ ych 3D vazb´ach, podrobnˇeji pops´ano v pr´aci). Napˇr. pro kmit´an´ı v jednom smˇeru lze redukovat rozmˇer u ´lohy na polovinu.

16

Poˇ c´ıtaˇ cov´ a implementace modelu Sch´ema realizace modelu zn´azorˇ nuje obr. 1. V preprocessingov´e ˇca´sti je pro tvorbu geometrie 1. Geometry of resonator 2. Setting boundary conditions 3. Building mesh

GMSH manually GMSH

*.msh file boundary_conditions.dta file

1. Computation of global matrices 2. Introduction of boundary conditions 3. Solving the eigenvalue problem freq.dta file coupling_coefficient.dta file

c++ code Arpack & c++ code

*.pos file

1. Selection of dominant modes 2. Vizualization of selected modes

c++ code GMSH

Obr´azek 1: Sch´ema implementace modelu a jeho ˇc´asti. a s´ıtˇe pouˇzit volnˇe dostupn´ y software GMSH [33]; ten tak´e slouˇz´ı k vizualizaci v´ ysledk˚ u. Pro sestaven´ı matic, zad´an´ı okrajov´ ych podm´ınek, identifikaci m´od˚ u kmit˚ u a selekci dominantn´ıch m´od˚ u byl vytvoˇren p˚ uvodn´ı programov´ y modul v jazyce C++. Pro ˇreˇsen´ı zobecnˇen´eho probl´emu vlastn´ıch ˇc´ısel je pouˇzit IRA v implementaci z volnˇe dostupn´e knihovny ARPACK [28] v jazyce FORTRAN. Informace mezi jednotliv´ ymi moduly jsou pˇred´av´any pomoc´ı textov´ ych soubor˚ u (podrobnˇejˇs´ı popis implementace naleznete v pr´aci).

17

Testov´ an´ı modelu Kmit´ an´ı kˇ remenn´ e tyˇ cinky Testov´an´ı kompletn´ıho modelu bylo provedeno na u ´loze kmit´an´ı kˇremenn´e tyˇcinky. Byly zvoleny pod´eln´e a torzn´ı m´ody kmit´an´ı a v´ ysledky byly porovn´any s analytick´ ym ˇreˇsen´ım. Geometrie a okrajov´ e podm´ınky. Geometrie, diskretizace a okrajov´e podm´ınky pro testovac´ı u ´lohu jsou zobrazeny na obr´azc´ıch n´ıˇze. Rozmˇery rezon´atoru byly zvoleny l = 0.01 m, a, b = 0.001 m.

+

b u=0

electrodes

-

u=0

a l

Obr´azek 2: Piezoelektrick´a tyˇcinka s um´ıstˇen´ım elektrod generuj´ıc´ım pod´eln´e a torzn´ı kmity a pˇredepsan´e okrajov´e podm´ınky.

Obr´azek 3: Pˇr´ıklad diskretizace rezon´atoru. V´ ysledky modelu dokumentuj´ı dobrou shodu s analytick´ ym ˇreˇsen´ım.

Pod´ eln´ e kmity Pro analytick´e vyj´adˇren´ı h-t´e rezonanˇcn´ı frekvence jsme pouˇzili vzorec [26] r h 1 , h = 1, 2, 3, ... fh = 2l %s11

18

Tabulka ukazuje porovn´an´ı analytick´eho a numerick´eho ˇreˇsen´ı pro prvn´ıch 6 frekvenc´ı. analyticky spoˇcten´e frekvence numericky spoˇcten´e frekvence relativn´ı odchylka(%) f1 = 271803 Hz f2 = 543606 Hz f3 = 815409 Hz f4 = 1087212 Hz f5 = 1359015 Hz f6 = 1630818 Hz

286205 Hz 572561 Hz 859212 Hz 1146292 Hz 1433952 Hz 1722252 Hz

+5.3 +5.3 +5.4 +5.4 +5.5 +5.6

Vizualizace v´ ysledn´ ych kmit˚ u je zn´azornˇena na obr´azku 4. Parametry v´ ypoˇctu byly2 : S´ıt’: poˇcet uzl˚ u = 4258, poˇcet prvk˚ u = 21326 Matice: A 7476 × 7476, 102446 nenulov´ ych prvk˚ u v horn´ı troj´ uheln´ıkov´e ˇca´sti, B 21326 nenulov´ ych prvk˚ u v horn´ı troj´ uheln´ıkov´e ˇca´sti ym ˇreˇsen´ım. Relativn´ı odKoment´ aˇ r v´ ysledk˚ u. Model prok´azal dobrou shodu s analytick´ chylky kolem 5% mohou b´ yt zp˚ usobeny v´ıce faktory - model zahrnuje vazby v rezon´atoru ve vˇsech smˇerech, zat´ımco analytick´e ˇreˇsen´ı ne (uvaˇzuje se jen koeficient s11 ). D´ale je model citliv´ y na zmˇenu vstupn´ıch materi´alov´ ych parametr˚ u. Pˇri zmˇenˇe materi´alov´ ych tenzor˚ u v ˇr´adu procent se tato zmˇena ˇr´adovˇe prom´ıtne i do v´ ysledk˚ u. Pouˇzit´e materi´alov´e tenzory jsou uvedeny v pr´aci.

Torzn´ı kmity Pro analytick´e vyj´adˇren´ı h-t´e rezonanˇcn´ı frekvence torzn´ıch kmit˚ u byl pouˇzit vzorec [26] s r 2 ab h G a ³a´ p fh = 1 − χ , h = 1, 2, 3, ... 2l % 1 + ( ab )2 b b ¡ ¢ kde G = 40973043485.312064 je efektivn´ı torzn´ı modul a hodnota funkce χ ab = 0.5777787. Srovn´an´ı analytick´eho a numerick´eho ˇreˇsen´ı je obsaˇzeno v n´asleduj´ıc´ı tabulce. analyticky spoˇcten´e frekvence numericky spoˇcten´e frekvence relativn´ı odchylka(%) f1 = 180668 Hz f2 = 361336 Hz f3 = 542004 Hz f4 = 722672 Hz f5 = 903340 Hz f6 = 1084008 Hz f7 = 1264676 Hz f8 = 1445344 Hz f9 = 1626012 Hz

171615 Hz 356761 Hz 532879 Hz 707107 Hz 853750 Hz 1024242 Hz 1189933 Hz 1361634 Hz 1627816 Hz

−5 −1.3 −1.7 −2.2 −5.5 −5.5 −6 −5.8 +0.1

Vizualizace kmit˚ u je zn´azornˇena na obr´azku 5. Parametry v´ ypoˇctu byly3 : S´ıt’: poˇcet uzl˚ u = 4258, poˇcet prvk˚ u = 21326 Matice: A 8606 × 8606, 223087 nenulov´ ych prvk˚ u v horn´ı troj´ uheln´ıkov´e ˇca´sti, B 49026 nenulov´ ych prvk˚ u v horn´ı troj´ uheln´ıkov´e ˇca´sti 2 3

Byla pouˇzita redukce na smˇer osy x. Byl ˇreˇsen kompletn´ı 3D probl´em.

19

Koment´ aˇ r v´ ysledk˚ u. Relativn´ı odchylky mohou b´ yt zp˚ usobeny faktory popsan´ ymi v´ yˇse. Nav´ıc algebraick´a u ´loha byla rozmˇernˇejˇs´ı neˇz v prvn´ım pˇr´ıpadˇe, coˇz s sebou m˚ uˇze n´est vˇetˇs´ı numerickou chybu. Nicm´enˇe, pr˚ umˇern´a relativn´ı odchylka je menˇs´ı neˇz v u ´loze pod´eln´ ych kmit˚ u (tento fakt lze pˇriˇc´ıst i tomu, ˇze analytick´ y vzorec jiˇz zahrnuje 2D vazbu - koeficienty s55 a s66 obsaˇzen´e v modulu G).

20

Obr´azek 4: Prvn´ıch 6 pod´eln´ ych m´od˚ u kmit˚ u piezoelektrick´e tyˇcinky. 21

Obr´azek 5: Prvn´ıch 8 torzn´ıch m´od˚ u kmit˚ u piezoelektrick´e tyˇcinky.

22

Aplikace modelu na re´ aln´ eu ´ loze – kmit´ an´ı planparaleln´ıho kˇ remenn´ eho rezon´ atoru Model byl aplikov´an na u ´loze nalezen´ı dominantn´ıch m´od˚ u kmit˚ u tlouˇst’kovˇe stˇriˇznˇe kmitaj´ıc´ıho planparaleln´ıho kˇremenn´eho rezon´atoru ve smˇeru osy x ˇrezu krystalu. Tento rezon´ator se pr˚ umyslovˇe vyr´ab´ı v podniku Krystaly, a.s., Hradec Kr´alov´e. Geometrie a ilustraˇcn´ı diskreti-

R mounting

u=0

rr

mounting u=0

electrodes

h

Obr´azek 6: Popis planparaleln´ıho rezon´atoru a ilustrativn´ı s´ıt’. zace rezon´atoru jsou zn´azornˇeny na obr. 6. Jedn´a se o kruhovou destiˇcku o polomˇeru R = 7 mm se dvˇema kruhov´ ymi elektrodami o polomˇeru r = 3.5 mm um´ıstˇen´ ymi ve stˇrdu horn´ı a spodn´ı plochy rezon´atoru. Rezon´ator je upevnˇen na protilehl´ ych stran´ach ve smˇeru osy x. Tlouˇst’ka rezon´atoru je h = 0.3355 mm. Rezon´ator je vyroben z kˇremenn´eho krystalu ˇrezu AT 35.25◦ . Materi´alov´e tenzory jsou uvedeny v pr´aci. Parametry v´ ypoˇctu byly4 : S´ıt’: poˇcet uzl˚ u = 7360, poˇcet prvk˚ u = 31860 Matice: A 9771 × 9771, 100430 nenulov´ ych prvk˚ u v horn´ı troj´ uheln´ıkov´e ˇca´sti, B 47038 nenulov´ ych prvk˚ u v horn´ı troj´ uheln´ıkov´e ˇca´sti V´ ysledky experiment´aln´ıho mˇeˇren´ı jsou pˇrevzaty z v´ yvojov´eho oddˇelen´ı firmy Krystaly, a.s. V´ ystup z mˇeˇren´ı je zn´azornˇen na obr. 8) se zv´ yraznˇen´ ymi hodnotami rezonanˇcn´ıch frekvenc´ı.

V´ ysledky a diskuze N´asleduj´ıc´ı tabulka 1 obsahuje porovn´an´ı namˇeˇren´ ych a spoˇcten´ ych hodnot pro prvn´ı tˇri rezonanˇcn´ı frekvence hledan´ ych m´od˚ u kmit˚ u. Pr˚ umˇern´a odchylka numerick´ ych v´ ysledk˚ u od namˇeˇren´ ych hodnot ˇcin´ı asi 15%, ale relativn´ı odstupy jednotliv´ ych dominantn´ıch m´od˚ u jsou dobˇre 4

Byla pouˇzita redukce na smˇer osy x.

23

rezonanˇcn´ı frekvence 1. (f1 ) 2. 3.

namˇeˇren´a (kHz) 4962 5067.5 = 1.02 f1 5102.5 = 1.03 f1

spoˇcten´a (kHz) 4210.843 4246.481 = 1.01 f1 4331.652 = 1.03 f1

relativn´ı odchylka 15.1% 16.2% 14.7%

Tabulka 1: Srovn´an´ı namˇeˇren´ ych a spoˇcten´ ych dominantn´ıch rezonanˇcn´ıch frekvenc´ı tlouˇst’kovˇe stˇriˇzn´ ych kmit˚ u. zachov´any. Pˇri hled´an´ı d˚ uvod˚ u t´eto odchylky lze odk´azat na diskuze k v´ ysledk˚ um u testovac´ıch u ´loh - model je citliv´ y na vstupn´ı materi´alov´e a geometrick´e parametry, nav´ıc se jedn´a o hled´an´ı m´od˚ u kmit˚ u nepomˇernˇe sloˇzitˇejˇs´ıch neˇz u testovac´ıch u ´loh. Byly nalezeny spr´avn´e dominantn´ı m´ody kmit˚ u. Graf na obr. 9 ukazuje rozloˇzen´ı koeficient˚ u elektromechanick´e vazby pro prvn´ıch 400 spoˇcten´ ych rezonanˇcn´ıch frekvenc´ı se zv´ yraznˇen´ım tˇr´ı hledan´ ych dominantn´ıch frekvenc´ı. Obr´azek ukazuje rozloˇzen´ı amplitud u spoˇcten´ ych m´od˚ u kmit˚ u a jejich porovn´an´ı s teoretick´ ym stavem. Pˇres jist´e nepˇresnosti, zˇrejm´e hlavnˇe u tˇret´ıho dominantn´ıho m´odu, lze jednotliv´e m´ody jednoznaˇcnˇe identifikovat s jejich teoretick´ ymi protˇejˇsky.

140

5400

120

5200

100

5000

80

4800

60

4600

40

4400

20

4200

0 0,2684

0,2768

0,2852

0,2936

0,3020

0,3103

0,3187

0,3271

frequency (kHz)

frequency separation (kHz)

ˇuje sledovat Moˇ zn´ e pouˇ zit´ı v´ ysledk˚ u pˇ ri procesu n´ avrhu rezon´ ator˚ u. Model umoˇzn chov´an´ı rezon´atoru v z´avislosti na zmˇen´ach tvaru ˇci um´ıstˇen´ı elektrod. Jako pˇr´ıklad, graf na obr. 7 ukazuje z´avislost dominant´ıch ´rezonanˇcn´ıch frekvenc´ı a jejich vz´ajemn´eho odstupu na zmˇenˇe tlouˇst’ky rezon´atoru. Model dobˇre zachycuje line´arn´ı z´avislost rezonanˇcn´ıch frekvenc´ı na tlouˇst’ce. Podobn´e v´ ysledky mohou b´ yt pouˇzity v procesu n´avrhu rezon´ator˚ u, napˇr. v optimalizaˇcn´ım procesu, kdy hled´ame takov´ y tvar a rozmˇery rezon´atoru, pro kter´e budou maximalizovan´e odstupy jednotliv´ ych dominantn´ıch frekvenc´ı ˇci odstupy dominantn´ıch frekvenc´ı od ostatn´ıch frekvenc´ı.

f_1 f_2 f_3 f_2-f_1 f_3-f_2 f_3-f_1

4000 0,3355

thickness (mm)

Obr´azek 7: Graf z´avislosti dominantn´ıch rezonanˇcn´ıch frekvenc´ı a jejich odstup˚ u na tlouˇst’ce rezon´atoru.

24

4 927 500

5 227 500

5 077 500 30 000

Obr´azek 8: Dominantn´ı rezonanˇcn´ı frekvence nalezen´e z mˇeˇren´ı.

120

electromechanical coupling coefficient

100

80

60

40

20

0 0

500

1000

1500

2000

2500

3000

3500

4000

4500

5000

5500

frequency (kHz)

Obr´azek 9: Graf koeficient˚ u elektromechanick´e vazby se zv´ yraznˇen´ım dominantn´ıch frekvenc´ı.

25

Obr´azek 10: Spoˇcten´e dominantn´ı m´ody a srovn´an´ı s teoretick´ ym rozloˇzen´ım amplitud. 26

Z´ avˇ er Dizertaˇcn´ı pr´ace pˇrin´aˇs´ı v´ ysledky modelov´an´ı rezonanˇcn´ıch charakteristik piezoelektrick´ ych rezon´ator˚ u. Diskretizace fyzik´aln´ıho modelu, zaloˇzen´a na metodˇe koneˇcn´ ych prvk˚ u, vede na zobecnˇen´ y probl´em vlastn´ıch ˇc´ısel, jehoˇz ˇreˇsen´ım lze z´ıskat rezonanˇcn´ı frekvence piezoelektrick´eho rezon´atoru. K ˇreˇsen´ı algebraick´e u ´lohy je pouˇzito Krylovovsk´ ych metod, jmenovitˇe implicitnˇe restartovan´eho Arnoldiho algoritmu. Poˇc´ıtaˇcov´a implementace modelu je pouˇziteln´a pro urˇcov´an´ı dominantn´ıch rezonanˇcn´ıch frekvenc´ı rezon´ator˚ u r˚ uzn´ ych tvar˚ u. Testovac´ı u ´lohy prok´azaly pouˇzitelnost modelu, pˇriˇcemˇz vypoˇcten´e v´ ysledky odpov´ıdaly analytick´emu ˇreˇsen´ı. Model, aplikovan´ y na re´alnou u ´lohu tlouˇst’kovˇe-stˇriˇzn´ ych kmit˚ u planparaleln´ıho kˇremenn´eho rezon´atoru, pˇrin´aˇs´ı kvalitativn´ı shodu s mˇeˇren´ım nalezen´ım odpov´ıdaj´ıc´ıch dominantn´ıch m´od˚ u pˇri zachov´an´ı odpov´ıdaj´ıc´ıch relativn´ıch odstup˚ u jednotliv´ ych frekvenc´ı. Odchylka od mˇeˇren´ı (asi 15%) ukazuje na omezen´ı modelu, ovˇsem po kalibraci na danou u ´lohu lze model prakticky pouˇz´ıvat. Citlivost modelu vyˇzaduje pˇresn´e zad´an´ı vˇsech vstupn´ıch parametr˚ u (zejm´ena materi´alov´ ych tenzor˚ uodchylka v zad´an´ı se ˇr´adovˇe prom´ıt´a do v´ ysledk˚ u). Model vykazuje spr´avnou odezvu na zmˇenu vstupn´ıch parametr˚ u (napˇr. zachycuje line´arn´ı z´avislost rezonanˇcn´ı frekvence na tlouˇst’ce rezon´atoru). Je zde nˇekolik oblast´ı pro budouc´ı pr´aci. Hlavn´ı je vytvoˇren´ı optimalizaˇcn´ıho modulu, jenˇz by byl pouˇziteln´ y pˇri procesu n´avrhu rezon´ator˚ u s poˇzadovan´ ymi vlastnostmi. To bude vyˇzadovat nˇekoliker´e opakov´an´ı ˇreˇsen´ı algebraick´e u ´lohy. S touto u ´lohou m˚ uˇze vyvstat poˇzadavek na rychlejˇs´ı poˇc´ıtaˇcovou implementaci, napˇr. pomoc´ı paraleln´ı implementace. Jak bylo zm´ınˇeno v pr´aci, proti rychlosti v´ ypoˇctu stoj´ı tak´e potˇreba rozmˇern´ ych u ´loh, pokud chceme vyj´adˇrit sloˇzit´e m´ody kmit˚ u. Moˇznou cestou k ˇreˇsen´ı rozmˇernˇejˇs´ıch u ´loh, neˇz dosud, je pouˇzit´ı nˇekter´e z implementac´ı v´ıce´ urovˇ nov´ ych metod pro ˇreˇsen´ı u ´lohy vlastn´ıch ˇc´ısel (viz napˇr. [2], [5] nebo [7]).

27

Literatura [1] Allik H., Hughes T. J. R., Finite element method for piezoelectric vibration, International journal for numerical methods in engineering, Vol. 2, pp. 151-157, 1970. [2] Bennighof J.K., Lehoucq R.B., An automated multilevel substructuring method for eigenspace computation in linear elastodynamics, SIAM, Vol. 25, 2003. [3] Brenner S.C., Scott L.R., The mathematical theory of finite element methods Springer-Verlag NY, 1994. [4] Brepta R., Pust L., Turek F., Mechanick´e kmit´an´ı, Technical Guide, Vol. 71, Sobot´ales, Prague 1994. [5] Chao Yang, Solving large-scale eigenvalue problems in SciDAC applications, Technical report, Lawrence Berkeley National Laboratory, 2005. [6] Dulmet B., Current trends of FEM applied to piezoelectric resonators, ENSMM, L’institut FEMTO-ST, Besancon. [7] Fialko S.Y., An aggregation multilevel iterative solver with shift acceleration for eigenvalue analysis of large-scale structures, CMM-2003 – Computer Methods in Mechanics, Gliwice, 2003. [8] Hwang W.S., Park H.Ch., Closed form expression for higher-order electroelastic tetrahedral elements, AIAA Journal, Vol. 31, No. 5, 1993. [9] Lehoucq R., Maschhoff K., Sorensen D., Yang Ch.: ARPACK USERS GUIDE, Solution of large scale eigenvalue problems by implicitly restarted Arnoldi methods, 1997. [10] Lerch R., Simulation of piezoelectric devices by two- and three-dimensional finite elements, IEEE Transactions on Ultrasonics, Ferroelectric and Frequency Control, Vol. 37, No. 2, 1990. [11] M´ıka S., Kufner A., Parci´ aln´ı diferenci´ aln´ı rovnice I, SNTL Praha 1983. [12] Milsom R. F., Elliot D. T., Terry Wood S., Redwood M., Analysis and design of couple mode miniature bar resonator and monolothic filters, IEEE Transactions on Ultrasonics, Ferroelectric and Frequency Control, Vol. 30, pp. 140- 155, 1983. [13] Moetakef M. A., Lawrence K.L., Joshi S.P., Shiakolas S., Closed form expression for higher-order electroelastic tetrahedral elements, AIAA Journal, Vol. 33, No. 1, 1995. ˇas J., Hlava ´c ˇek I., Mathematical theory of elastic and elasto-plastic bodies: An [14] Nec introduction, 1. ed., Elsevier, Amsterdam 1981.

28

[15] Piefort V., Finite element modelling of piezoelectric active structures, Ph.D. thesis, FAS ULB, 2001. [16] Rektorys K., Variaˇcn´ı metody, Academia Praha 1989. [17] Samanta B., Ray M.C., Bhattacharyya R., Finite element model for active control of intelligent structures, AIAA Journal, Vol. 34, No. 9, 1996. ˇenko V.K., Izbrannye glavy teoretiˇceskoj fiziki, Voennogo izdatelstvo Ministerstva [18] Semenc oborony SSSR, Leningrad 1960. [19] Sorensen D., Implicitly restarted Arnoldi/Lanczos methods for large scale eigenvalue calculations, 1995. [20] Stewart G. W., Matrix algorithms, University of Maryland, College Park 1999. ˇ [21] Stoll I., Tolar J., Teoretick´ a fyzika, Czech Technical University, Prague 1981. [22] Tiersten H.F., Linear piezoelectric plate vibrations: Elements of the linear theory of piezoelectricity and the vibrations of piezoelectric plates, Plenum Press, New York 1969. [23] Trefethen L. N., Bau D., Numerical linear algebra, Siam, Philadelphia 1997. [24] Tzou H. S., Tseng C. I., Distributed modal identification and vibration control of continua: piezoelectric finite element formulation and analysis, Journal of Dynamic Systems, Measurement and Control, Vol. 113, pp. 500- 505, 1991. [25] Wang J., Yong Y.K., Imai T., Finite element analysis of the piezoelectric vibrations of quartz plate resonators with higher-order plate theory, Proceedings of 1997 IEEE International frequency control symposium, Orlando 1997. ´ J., Zelenka J., Pod´elnˇe a ploˇsnˇe stˇriˇznˇe kmitaj´ıc´ıpiezoelektrick´e rezon´ [26] Tichy atory ze ˇ Cas. ˇ synthetick´eho kˇremene, Cs. Fys. 10, pp. 328 - 332, 1960. [27] Zelenka J., Piezoelektrick´e rezon´ atory a jejich pouˇzit´ı v praxi, Academia, Praha 1981. [28] Lehoucq R., Maschhoff K., Sorensen D., Yang Ch., www.caam.rice.edu/software/ARPACK/ [29] Paige C. C., Saunders M. A., http://www.stanford.edu/group/SOL/software/symmlq.html [30] Marques O., SKYPACK User’s Guide, http://crd.lbl.gov/∼osni/#Software [31] Marques O., http://crd.lbl.gov/∼osni/#Software [32] http://www.netlib.org/lapack [33] http://www.geuz.org/gmsh/ [34] http://www.geuz.org/gmsh/doc/texinfo/gmsh.htm [35] http://www.krystaly.cz

29

[36] http://www.sensortech.ca Vybran´ e publikace k t´ ematu dizertace, napsan´ eˇ ci participovan´ e autorem: ´ lek P., Modelov´ [37] Ra an´ı piezoelektrick´ych jev˚ u, Diploma thesis, FJFI, Czech Technical University, Prague 2001. ´ k J., Ra ´ lek P., Application of FEM in modelling of the resonance [38] Maryˇ ska J., Nova characteristics of piezoelectric resonators, Proceedings of Algoritmy 2002, Slovakia, pp. 215-222, 2002. ´ k J., Ra ´ lek P., Finite element model of piezoelectric resonator, [39] Maryˇ ska J., Nova Current Trends in Scientific Computing, Contemporary Mathematics, Vol. 329, pp. 263270, 2003. ´k J., Ra ´lek P., Modelling of resonance characteristic of piezoelectric [40] Maryˇ ska J., Nova resonators - experimental experience, Proc. 6th ECMS 2003, Liberec, pp.123-128, 2003. ´k J., Ra ´ lek P., FEM model of piezoelectric continuum and selected [41] Maryˇ ska J., Nova applications, Proceedings of the IMAMM 2003, Roˇznov p. R., 2004. ´ lek P., Modelling of Piezoelectric Materials, Proceedings of the IX. PhD. Conference [42] Ra ˇ 2004. ICS, AV CR, ´ lek P., Generalized eigenvalue problem in modelling of the piezoelectric resonators, [43] Ra Proceedings of Seminar on numerical analysis SNA ‘05, Ostrava, 2005. ´ k J., Ra ´ lek P., FEM modelling of the resonance frequencies of the [44] Maryˇ ska J., Nova planparallel quartz resonator, Proceedings of ECMS 2005, Toulouse, 2005.

30