MIKROMECHANIKAI ÉRZÉKELŐK ELEKTROMOS-, MECHANIKAI- ÉS TERMIKUS MODELLEZÉSE
Doktori értekezés
Vízváry Zsolt okleveles gépészmérnök
Tudományos Vezető: Dr. Kovács Ádám
Budapest, 2005. április
i
BUDAPESTI MŰSZAKI ÉS GAZDASÁGTUDOMÁNYI EGYETEM GÉPÉSZMÉRNÖKI KAR Szerző neve:
Vízváry Zsolt
Értekezés címe: Mikromechanikai érzékelők elektromos-, mechanikai- és termikus modellezése Témavezető neve (ha volt): Dr. Kovács Ádám Értekezés benyújtásának helye (Tanszék, Intézet): Budapesti Műszaki és Gazdaságtudományi Egyetem, Műszaki Mechanikai Tanszék Dátum: 2005. április 5. Bírálók:
Javaslat:
nyilvános vitára igen/nem bíráló neve nyilvános vitára igen/nem bíráló neve nyilvános vitára igen/nem bíráló neve (ha van)
A bíráló bizottság javaslata: Dátum:
(név, aláírás) a bíráló bizottság elnöke
ii
Nyilatkozat
Alulírott Vízváry Zsolt kijelentem, hogy ezt a doktori értekezést magam készítettem és abban csak a megadott forrásokat használtam fel. Minden olyan részt, amelyet szó szerint, vagy azonos tartalomban, de átfogalmazva más forrásból átvettem, egyértelműen, a forrás megadásával megjelöltem. A dolgozat bírálatai és a védésről készült jegyzőkönyv a Budapesti Műszaki és Gazdaságtudományi Egyetem Gépészmérnöki Karának dékáni hivatalában elérhető.
Budapest, 2005. április 5.
Aláírás
iii
Electrical-, Mechanical- and Thermal modelling of Micromechanical Structures
Abstract The topic of this work is modelling coupled mechanical problems, which occurs in MEMS (Micro-Electro-Mechanical Systems). In the most cases one has to deal with mechanical, thermal and electric effects. In the first parts in my work the derivation and finite element discretization of piezothermoelastic equations are presented. I derive the coupled piezothermoelastic equations using a modified hyperbolic equation of heat conduction which is more adequate in micro scale. The equations are valid for homogenous, linear elastic, anisotropic material. Numerical calculations are carried out on modelling examples with an own finite element code. The next part shows the thermomechanical analysis of real micromechanical structures. The effects of geometrical parameters, material properties and boundary conditions are presented. Finally, elasticity modulus measurement of silicon-nitride is presented.
1
Tartalom
TARTALOM
1. Bevezetés....................................................................................................................................................... 2 2. Piezo-termo-elasztikus alapegyenletek.......................................................................................................... 8 2.1 Lineáris kapcsolt piezo-termo-elasztikus egyenletek .............................................................................. 8 2.2 Az általánosított lineáris piezo-termo-elasztikus egyenletek................................................................. 17 2.3 Kezdeti és peremfeltételek..................................................................................................................... 20 3. Piezo-termo-elasztikus problémák végeselemes megfogalmazása.............................................................. 23 3.1 Végeselemes egyenletek hely szerinti diszkretizációja ......................................................................... 23 3.1.1 Hely szerinti diszkretizáció a klasszikus hővezetési egyenlet esetén ............................................. 25 3.1.2 Az együttható mátrixok elemei a klasszikus hővezetési egyenlet esetén........................................ 32 3.1.3 Hely szerinti diszkretizáció a módosított hővezetési egyenlet esetén............................................. 34 3.1.4 Az együttható mátrixok elemei a módosított hővezetési egyenlet esetén....................................... 36 3.2 Végeselemes egyenletek idő szerinti diszkretizációja ........................................................................... 37 4. Végeselem program piezo-termo-elasztikus feladatokhoz (PiTE) .............................................................. 39 4.1 Elemtípus............................................................................................................................................... 39 4.2 Teszt példák........................................................................................................................................... 40 4.2.1 Mechanikai analízis ........................................................................................................................ 41 4.2.2 Termomechanikai analízis .............................................................................................................. 46 4.2.3 Piezoelektromos analízis ................................................................................................................ 50 4.3 Mintapélda............................................................................................................................................. 54 4.3.1 Mindkét végén befogott tartó.......................................................................................................... 54 5. Mikromechanikai szerkezetek kapcsolt modellezése .................................................................................. 61 5.1. Egy dimenziós kapacitív gyorsulásérzékelő ......................................................................................... 61 5.1.1. Mechanikai analízis ....................................................................................................................... 62 5.1.2. Termomechanikai analízis ............................................................................................................. 64 5.1.3. A mechanikai és termomechanikai eredmények összehasonlítása ................................................ 64 5.2. Gázérzékelő .......................................................................................................................................... 71 5.2.1. Négyhidas gázérzékelő szerkezet .................................................................................................. 72 5.2.2. Kéthidas gázérzékelő szerkezet ..................................................................................................... 77 6. Szilícium-nitrid rugalmassági moduluszának meghatározása ..................................................................... 81 6.1.Mérés..................................................................................................................................................... 81 6.1.1.Mérési elrendezés ........................................................................................................................... 81 6.1.2. A minták ........................................................................................................................................ 83 6.2. Számítások............................................................................................................................................ 84 6.2.1. Analitikus modell .......................................................................................................................... 84 6.2.2. Végeselem modell ......................................................................................................................... 85 6.2.3. Eredmények................................................................................................................................... 85 Összefoglalás................................................................................................................................................... 88 Summary ......................................................................................................................................................... 89 Tézisek ............................................................................................................................................................ 90 Irodalom .......................................................................................................................................................... 91 Publikációk...................................................................................................................................................... 98 Köszönetnyilvánítás ...................................................................................................................................... 101
1. Bevezetés
2
1. BEVEZETÉS A korszerű, nagy sebességű számítógépekkel ma már rendkívül sokféle folyamatirányítási és szabályozási feladat oldható meg. Ezeknek a feladatoknak az eredményes megoldásához feltétlenül szükséges, hogy az ellenőrizendő és szabályozni kívánt folyamatokat jellemző fizikai, kémiai paramétereket folyamatosan érzékeljük, mérjük. Ezekhez a mérésekhez szükség van érzékelőkre (szenzorok). A szabályozási paramétereket nemcsak mérni kell, de esetenként be is kell avatkozni a folyamatba. Ezeket az eszközöket, amelyeket fordított működésű érzékelőknek is tekinthetünk, beavatkozóknak (aktuátorok) nevezzük. Az érzékelők és beavatkozók kutatása, fejlesztése napjainkra kulcsfontosságúvá vált. A félvezető technika rohamos fejlődése lehetővé tette az egyre nagyobb mértékű miniatürizálást. Ez az érzékelők és beavatkozók szempontjából rendkívül fontos, hiszen minél kisebb a mérőeszköz, annál kevésbé zavarja meg a szabályozni kívánt folyamatot. A félvezető technika esetében pedig még arra is lehetőség van, hogy a kiértékelő áramkör is elkészíthető egy chipen belül. Az érzékelők és beavatkozók többféle működési elven alapulnak: pl. mechanikai, optikai, mágneses, stb. Ezeken a fő kategóriákon belül azonban újabb működési elvek különböztethetőek meg: pl. a mechanikai érzékelőkön belül beszélhetünk piezorezisztív, kapacitív, piezoelektromos, rezgőelemes stb. érzékelőkről. Dolgozatomban kizárólag a mechanikai érzékelőkkel foglalkozom, azon belül is elsősorban a piezorezisztív, piezoelektromos és kapacitív működési elvűekkel. A piezorezisztív és kapacitív érzékelők elsődleges alapanyaga a szilícium. A szilícium rendkívül jó mechanikai tulajdonságokkal rendelkező anyag. A szilícium egykristály szinte tökéletesen lineárisan rugalmasan, rideg anyagként viselkedik szobahőmérsékleten, folyás nélkül törik egészen 600 oC-ig. Magasabb hőmérsékleten a törés előtt kis mértékű folyás tapasztalható [1]. Emellett elektromos tulajdonságai is jók, ezért is kedvelt alapanyag a mikromechanikában. Szilícium alapú érzékelők lehetnek például nyomásmérők, erőmérő cellák, gyorsulásérzékelők, amelyekkel a hétköznapok során is találkozhatunk. Az ilyen eszközöket a nemzetközi irodalomban MEMS-nek nevezik (Micro Electro-Mechanical System). Az elektromechanikai érzékelők működése azon alapul, hogy a különböző energiákat átalakítják – mechanikai energiából elektromos energiát, illetve fordítva -, ezért fejlesztésük kutatása több tudomány - anyagtudomány, elektromosságtan, hőtan, mechanika, fizika, kémia - határterületére esik. Kapacitív érzékelők esetében a mechanikai jellegű információ az elmozdulás, amit elektromos jellé kívánunk alakítani. Ezt úgy tehetjük meg, hogy készítünk egy kondenzátort, melynek egyik fegyverzete a mozgó elem, ami általában membrán, de lehet bármilyen más felfüggesztett struktúra, a másik pedig egy rögzített felület, ami a mozgó felülettel szemben van. Így az elmozdulás változás a kapacitás megváltozását okozza. Minél nagyobb az elmozdulás, annál nagyobb a kapacitás változása is. Számítási nehézséget okozhat az eszközök nagy részénél az, hogy a mozgó rész rugalmas, ami miatt
1. Bevezetés
3
az elmozdulás nem konstans a lemez pontjaiban, így a kapacitás az elmozdulásnak nem lineáris függvénye. Gondoljunk például a membránokra, ahol a membrán közepén a legnagyobb a lehajlás. Ez a tervezhetőség szempontjából is bonyolítja a helyzetet, hiszen nem síkkondenzátorról van szó, így a várható kapacitásváltozás kiszámításához is numerikus szimulációra lehet szükség. Kapacitív típusú érzékelőket beavatkozóként is lehet alkalmazni. Az elektromos feszültség hatására kialakuló térerősség a mozgó fegyverzet elmozdulását okozza. Piezorezisztív érzékelők esetén olyan érzékelőkről beszélünk, ahol ellenállást (piezorezisztor) alakítanak ki az érzékelő azon helyén, ahol a szerkezetben mechanikai feszültségek keletkezhetnek. A piezoellenállás működése a nyúlásmérő bélyegekéhez hasonló. Mechanikai feszültség változás hatására ellenállás változás következik be. Az ellenállásokat általában Wheatstone-hídba szokták kapcsolni. A piezorezisztivitás nem használható beavatkozókban. Piezoelektromos jelenségről akkor beszélünk, ha mechanikai feszültség hatására, elektromos feszültség, illetve áram keletkezik a szerkezetben, vagy fordítva. Ebből adódóan a piezoelektromos eszközök beavatkozóként is működhetnek. Problémát okozhat viszont az, hogy hosszabb idejű statikus terhelés esetén felléphet szökési áram, azaz az elektromos jel időben változhat annak ellenére, hogy a mechanikai terhelés változatlan. Kezdetben az elektromechanikai érzékelők fejlesztése – talán amiatt, hogy félvezető technikai alapokon nyugszik – a fizikusok és villamosmérnökök szakterülete volt. Ennek eredményeként ezeket az érzékelőket – annak ellenére, hogy valamilyen mechanikai jellemző mérésére fejlesztették – mechanikai szempontból nem vizsgálták. Ennek még ma is érezhető az a kedvezőtlen hatása, hogy az érzékelők és beavatkozók anyagjellemzőire vonatkozóan korlátozottak, olykor ellentmondásosak az irodalomban található adatok, noha mára már számos publikációt találhatunk, amely a mikromechanikai szerkezetek mechanikai vizsgálatával foglalkozik. Az érzékelők és beavatkozók vizsgálatát nehezíti az is, hogy gyakran nem csupán egy mechanikai problémát kell vizsgálni, hanem valamilyen kapcsolt feladatot is. Hogy pontosan milyet, az az adott eszköztől függ. Attól függően, hogy mely mezők vannak kapcsolva, beszélhetünk elektromechanikai (piezoelektromos), termomechanikai vagy mindhárom mező kapcsoltsága és rugalmas alakváltozás esetén ún. piezo-termo-elasztikus (piezohőrugalmas) feladatról. A piezo-termo-elasztikus konstitutív egyenletek levezetését először MINDLIN-nél [2] találhatjuk meg. Az egyenletek levezetése, összefoglalása NOWINSKI [3] könyvében is szerepel. Ő ezekre az egyenletekre a piezo-termo-elasztikus elnevezést használja. Az elnevezések tekintetében az irodalom még nem egységes e tekintetben. A legelterjedtebb a piezo-termo-elasztikus elnevezés, de ismeretesek még a termo-elektro-mechanikai, termo-elektroelasztikus vagy termo-piezoelektromos elnevezések is. ALTAY és DÖKMECI [4] különbséget tesz a piezo-termo-elasztikus és a termo-piezo-elasztikus elnevezések között. Ők a piezo-termo-elasztikus egyenleteket a termo-piezo-elektromos egyenletek speciális esetének tekintik, ahol a hőmérséklet mező
1. Bevezetés
4
nincs a piezoelektromos egyenletekhez kapcsolva. Cikkükben tévesen hivatkoznak NOWINSKI-re is, mivel a [3] könyvben leírt egyenletek összekapcsolják a hőmérsékletmezőt az elektromos és mechanikai mezőkkel. Dolgozatomban a piezo-termoelasztikus, illetve piezo-hőrugalmas elnevezéseket fogom használni NOWINSKI nyomán. Piezo-termo-elasztikus feladatok anyagegyenletei már az 1970-es évektől ismertek, azonban ezekkel a feladatokkal az 1990-es évekig nem sokan foglalkoztak. A 90-es évektől azonban – talán épp a szenzortechnika fejlődése miatt – az ilyen publikációk száma rohamosan nő. Kezdetben főleg analitikus megoldási lehetőségekkel foglalkoznak (pl. végtelen hosszú szalag, félvégtelen tér), de az ezredfordulóhoz közeledve már található egy-két olyan cikk is, ami ezeknek az egyenleteknek a végeselemes megoldását tárgyalja. A kapcsolt piezo-termo-elasztikus feladatok analitikus megoldására speciális példák vannak az irodalomban. MINDLIN [5] nagyfrekvenciájú rezgéseket vizsgál kétdimenziós egyenletekkel piezoelektromos kristály lemezeken a termikus hatást is figyelembe véve. DUBE és társai [6] pontos analitikus megoldást ad egy végtelen hosszú, koncentrált nyomatékkal terhelt piezo-termo-elasztikus ortotróp szalagra. LEE és SARAVANOS [7] kompozit héjakat vizsgál a végeselem módszer segítségével, de cikkükben a hőmérséklet csak az elmozdulás és elektromos mezőkhöz tartozó terhelés vektorban jelenik meg, tehát valójában csak két mező kapcsolásáról van szó. ASHIDA és TAUCHERT [8] piezo-termoelasztikus hengeres lemezek síkfeszültségi problémáival foglalkozik. Hivatkozott cikkük korábbi cikkeiknél általánosabb eseteket vizsgál. STAM és CARMAN [9] koncentrikusan egymáshoz csatolt piezoelektromos hengerek termoelektromechanikai viselkedését vizsgálja. CHEN és SHEN [10] hengerszimmetrikus piezo-termo-elasztikus héjat vizsgál analitikus módszerekkel. QIN [11] repedéseket vizsgál piezoelektromos félsíkon a peremelem módszerrel. QIN [12], [13], [14] piezo-termo-elasztikus Green-függvény előállításával foglalkozik lemezekre különböző nyílásokkal. BELYAEV [15] a hővezetési egyenletet és dinamikus határérték problémát vizsgál piezo-termo-elasztikus anyagokban. TARN [16] és MELNIK [17] numerikus módszerekkel közelítik meg a problémát, de nem a végeselem módszerrel. Végeselemes megfogalmazást találhatunk FUNG és társai [18], valamint KOKO és társai [19] cikkében, de egyikben sem közlik az egyenletek alapjait. Részletes levezetést közöl BOURANTA és társai [20]. Ők virtuális elmozdulásokat és variációs elveket használnak a végeselemes egyenletek levezetéséhez. GÖRNANDT és GABBERT [21] szintén variációs elveket használnak az egyenletek felírásához. Az utóbbi két cikk gyakorlatilag ugyanazt a mátrix differenciálegyenletet eredményezi (BOURANTA-ék cikke [20] tartalmaz egy kis elírást), azonban eredményük eltér mind a [18], mind a [19] cikk végeselemes felírásától. Dolgozatomban nem variációs elveket fogok használni a végeselemes egyenletrendszer felírásához, hanem a Galjorkin súlyozott maradvány módszert. A tisztán mechanikai és termomechanikai feladatok megoldása elvben ma már nem okoz nagy problémát. A legtöbb végeselem programcsomag tartalmaz olyan modulokat, amelyek segítségével ezek a feladatok megoldhatóak. Természetesen ezen a területen is
1. Bevezetés
5
adódnak olyan feladat típusok, amelyeket ezek a programok nem fednek le (pl. általánosított hővezetési egyenlet). BALLA [22] kandidátusi értekezésében nagyon alaposan levezeti általános esetre a hőrugalmasságtan végeselemes egyenleteit. Dolgozatomban nagy mértékben támaszkodom az ő munkájára. A piezoelektromos esetekre is vannak már modulok, sőt az ANSYS [23] végeselem programnál már bizonyos hármas kapcsoltság is van, de ott valójában csak az elektromos és a termikus, illetve a termikus és a mechanikai mezők kapcsoltak páronként. FUNK és társai [24] háromdimenziós termo-elektro-mechanikai szimulációt ismertetnek MEMS eszközökre, de cikkükből kiderül, ők is csak páronként kapcsolták össze az egyenleteket. Annak ellenére, hogy a tisztán mechanikai és termomechanikai feladatok megoldhatóak a kereskedelmi végeselem programokkal, modellezés szempontjából a mikromechanikai eszközök modellezése tartogat még kihívásokat, hiszen szükségesek az anyagjellemzők. Sok esetben nem izotróp a feladat (pl. szilícium). CHEN és társai [25] cikkében azzal foglalkozik, hogy mi a hatása annak, hogyha a szilíciumot izotróp anyagnak tételezzük fel mikrostruktúrák tervezése során. Kérdéses az is, hogy azok az anyagjellemzők, amelyeket nagyobb próbatesteken meghatároznak, használhatóak-e kis méretekben, illetve vékony rétegek esetén. MAZZA és DUAL [26] mikrométer nagyságú, éles sarkokkal kialakított szilícium egykristály szerkezetek mechanikai viselkedését vizsgálja. PEARSON és társai [1] kisméretű szilícium kristályok deformációjával és törésével foglalkozik. WEI és társai [27] a mikroszerkezet hatását vizsgálják a szilícium vékonyrétegek hővezetése tekintetében. FANG és társai [28], [29] hőtágulási együtthatót határoznak meg mikroméretű vékonyrétegekben. Szintén a hőtágulási együttható meghatározásával foglalkozik ZIEBART, és társai [30]. A szerkezetek geometriája is fontos lehet. Általában különböző marási technológiákkal alakítják ki az eszközöket, de különböző rétegeket is képezhetnek a felületen. Ez különösen fontos a modellezés szempontjából, hiszen a különböző anyagú rétegek különböző anyagjellemzőik miatt okoznak problémát. Anyagjellemzők tekintetében az első publikáció, ami a szilíciummal, mint mechanikai anyaggal foglalkozik PETERSEN cikke [31]. Azonban a cikk nagy része nem az anyagjellemzők meghatározásával, hanem a szilícium megmunkálhatóságával foglalkozik. Jól szemlélteti, hogy félvezető technológiai módszerekkel milyen módon alakíthatóak ki membránok, bemélyedések, lyukak és hogyan befolyásolják a szilícium kristálytani jellemzői ezeket a technológiai folyamatokat. Az egyes technológiai lépések folyamán a szilíciumot szennyezik, doppolják. Ennek következtében az anyagjellemzői megváltozhatnak, illetve maradó feszültség keletkezik, azaz előfordulhat, hogy a kialakított szerkezet tehermentes állapotban sem feszültségmentes. Ezért szükséges lehet a szerkezetben e maradó feszültségek meghatározása. Ezzel foglalkozik PAN és HSU [32] mikrostruktúrák esetén. LEE és társai [33] a szennyező anyagok hatását tanulmányozza a Young-moduluszra. CHU, MEHREGANY [34] a maradó feszültségek vastagság menti eloszlását vizsgálja vékony
1. Bevezetés
6
p+-rétegekben1. DING, és társai [35] bórral adalékolt p+-rétegekben ébredő maradó feszültségekkel és a réteg mechanikai tulajdonságaival foglalkozik. A mikromechanikai technológiákkal előállított elektro-mechanikai érzékelők között a legelterjedtebbek a nyomásérzékelők és gyorsulásérzékelők. Még manapság is jelennek meg olyan cikkek ebben a témában, amelyek az érzékelő mechanikai viselkedésével, jellemzőivel mechanikai méretezésével nem foglalkoznak. Az érzékelő jellemzőit csak mérés útján állapítják meg. Ennek az lehet az oka, hogy a mikromechanikai eszközök előállítása nem jár nagy költségekkel egy már meglévő tiszta térrel rendelkező intézetben, üzemben, így néhány prototípus előállításával és teszteléssel próbálják megtalálni a megfelelő paramétereket. Nem tartalmaznak mechanikai számítást a következő cikkek [36] (1983), [37] (1989), [38] (1997), [39] (1998), [40], [41] (1999), [42], [43], [44], [45] (2000). Az élesedő piaci verseny azonban megköveteli a szerkezetek mechanikai tervezését, ellenőrzését is. Az ezredfordulóhoz közeledve egyre több az olyan publikáció, ami már tartalmaz valamiféle mechanikai számítást. Az első – általában gyorsulásérzékelővel kapcsolatos – publikációk analitikus képleteket tartalmaznak ez első sajátfrekvencia meghatározására. Ritkábban előfordul, hogy a lehajlást is meghatározzák analitikusan. ROYLANCE és ANGELL [46] egyik végén befogott tartó jellegű gyorsulásérzékelő előállításával foglalkozik. SEIDEL és CSEPREGI [47] is hasonló szerkezetek optimalizálást tűzi ki feladatául. Egy ennél bonyolultabb két tömeggel rendelkező gyorsulásérzékelő sajátfrekvenciáját vizsgálja SHENés társai [48]. ELGAMEL [49] nyomásérzékelő membránok lehajlás-kapacitásváltozás közötti kapcsolattal foglalkozik. A mikromechanikai szerkezetek modellezésére bonyolultabb geometriák esetén kézenfekvőnek látszik a végeselem módszer. CAO és társai [50] egy piezorezisztív membrán elkészítését és modellezését mutatja be. Ők is – és ez általában jellemző a mikromechanikai szerkezetekkel foglalkozó cikkekre – az ANSYS programot használták. Cikkükben a szerkezetben fellépő alakváltozásokat számították ki. CRAZZOLARA és társai [51] piezorezisztív gyorsulásérzékelő vizsgálatára használják a végeselem módszert, azonban a közölt eredmények csak kvalitatívak (csupán a feszültségkomponens előjelére vonatkozóan tesz megállapításokat). YOSHIDA és társai [52] egy kapacitív típusú gyorsulásérzékelő vizsgálatát végezte el. A számítások az elmozdulás-mezőre vonatkoznak, ugyanis ebből számolható ki a kapacitásváltozás. Szintén gyorsulásérzékelő modellezését találhatjuk meg LÜDTKE és társai [53] munkájában. Ez szintén az elmozdulás mező meghatározására terjed ki, az eredményeket analitikus számításokkal is összevetik. LI és társai [54] is az elmozdulásokat számítja ki egy giroszkóp modellezése során. CHOUAF és társai [55] már a mikromechanikai módszerekkel előállított membránok kritikus pontjaiban vizsgálja a feszültségeket. Többfajta végeselem modellt is készítenek a feladat modellezésére. XIAO és társai [56] gyorsulásérzékelő modális analízisét végzi el. PLAZA és társai [57] cikkében már egy átfogó mechanikai vizsgálatot láthatunk. Mind az elmozdulás mezőt, mind a feszültségmezőt meghatározták, ezen kívül a 1
Félvezető anyagok esetén p típusúra adalékolt réteg.
1. Bevezetés
7
szerkezet sajátfrekvenciáit is kiszámolták. BLASQUEZ és társai [58] pedig már egy nyomásérzékelő termomechanikai számításait mutatja be. Kapacitív típusú érzékelőknél a két fegyverzet között erőhatás léphet fel az elektromos térerő miatt. Ennek hatásával foglalkozik PUERS és LAPADATU [59]. LUCK és AGBA [60] piezoelektromos érzékelők és beavatkozók tervezési alapelveivel foglalkozik. A cikk tulajdonképpen a piezo-termo-elasztikus alapegyenleteket mutatja be. A kutatások ma már nemcsak a mikromechanikai eszközök tervezésével, hanem a megbízhatósági kérdéseivel is foglalkoznak. A szilícium törésmechanikai jellemzőivel foglalkozik HAUCH és társai [61]. SHARPE és társai [62] a mintadarab méretének hatását vizsgálja a Young-moduluszra és a törési szilárdságra. TAYEBI és TAYEBI [63] MEMSeszközök numerikus törésmechanikai analízisével foglalkozik. MUHLSTEIN, BROWN [64] és MESCHEDER, KRONAST és NAYCHUK [65] megbízhatósági és kifáradási jelenségeket vizsgálnak mikromechanikai eszközök esetében. Dolgozatom első részében áttekintést nyújtok a piezo-termo-elasztikus konstitutiv egyenletekről, azok levezetéséről. Galjorkin módszerével elvégzem a piezo-termoelasztikus egyenletek hely szerinti diszkretizációját mind a klasszikus, mind az általánosított hővezetési egyenlet felhasználásával. Elkészítettem egy végeselem programot, amely alkalmas piezo-termo-elasztikus statikai problémák megoldására. A programot FORTRAN nyelven írtam SMITH és GRIFFITHS [66] könyvének és mintaprogramjainak felhasználásával. Az általam kapott eredményeket összevetettem mind a könyvbeli programok, mind a COSMOS/M 2.0 [67] végeselem program eredményeivel. Szót ejtek még az időbeli diszkretizáció lehetőségeiről, hiszen a teljes mátrix differenciálegyenlet megoldásához erre is szükség van. Konkrét mikromechanikai szerkezetek mechanikai, illetve termomechanikai vizsgálatát végeztem el szintén a COSMOS/M 2.0, illetve az ANSYS programcsomaggal, mely az 5. fejezetben található.
8
2. Piezo-termo-elasztikus alapegyenletek
2. PIEZO-TERMO-ELASZTIKUS ALAPEGYENLETEK 2.1 Lineáris kapcsolt piezo-termo-elasztikus egyenletek A piezo-termo-elasztikus feladatok lineárisan rugalmas tulajdonságú anyagok esetén kerülnek elő, ezért az anyagtörvények levezetéséhez kihasználhatjuk ezt, és kis alakváltozásokat tételezünk fel. A piezo-termo-elasztikus egyenletek három – a mechanikai, a termikus és az elektromos – mező összekapcsolásával adódnak. Az anyagtörvények levezetéséhez szükségesek a három mezőt leíró egyenletek. Az elektromágneses feladatokat a Maxwell-egyenletek fogalmazzák meg. Ezek az egyenletek mind az elektromos, mind a mágneses viselkedést leírják. A piezo-termo-elasztikus egyenletekben azonban csak az elektromos viselkedésnek kell szerepelnie, így megtehetjük, hogy a mágneses tér hatását leegyszerűsítjük, vagy elhanyagoljuk. Az anyagtörvények levezetésének összefoglalását NOWINSKI [3] munkája alapján tehetjük meg. A levezetés során a szorzások jelei: Skaláris szorzat: ⋅ ( a ⋅ b ) Kétszeres skaláris szorzat: : ( A : B ) Vektoriális szorzat: × ( a × b ) Diadikus szorzat: o ( a o b ) Skalárral való szorzás esetén nem teszek szorzás jelet. ( ab ) A vektor és tenzor mennyiségeket vastag dőlt betűk jelölik a skalárokat normál álló betűk. A Maxwell-egyenletek:
∇× H = J +
∂D , ∂t (2.1)
∇× E = −
∂B , ∂t (2.2)
∇⋅B = 0,
(2.3) ∇⋅ D = γ , (2.4) ahol H – a mágneses térerősség, E – az elektromos térerősség, D – az elektromos eltolás vektor, B – a mágneses induktivitás,
9
2. Piezo-termo-elasztikus alapegyenletek
J – az áramsűrűség, γ – az elektromos töltés térbeli sűrűsége. A mechanikai tartalmat a mozgásegyenletek képviselik. A mozgásegyenlet felírásánál kihasználható az, hogy csak kis alakváltozásokat tételezünk fel, így nem teszünk különbséget a kezdeti és a pillanatnyi konfiguráció között. Minden az ismert kezdeti konfigurációban írható fel: && , σ ⋅ ∇ + F = ρu (2.5) σ =σ , T
(2.6) ahol σ - a Cauchy-féle feszültség tenzor, F – a térfogaton megoszló erőrendszer vektora, u – az elmozdulás vektor, ρ – sűrűség. Szintén a kis alakváltozás miatt a geometriai egyenlet a következő formában írható: ε=
1 (∇ o u + u o ∇ ) , 2 (2.7)
ahol ε – az alakváltozási tenzor. A termikus részt a hőáramra vonatkozó konstitutív összefüggések, illetve az entrópiára vonatkozó egyenletek felírásával építhetjük be az anyagtörvényekbe. Az anyagtörvényeknek ki kell elégíteniük a termodinamika I. és II. főtételét is. Ahhoz, hogy ez teljesüljön, célszerű a konstitutív egyenletek levezetésénél a főtételekből kiindulni. Az I. főtétel esetében a globális alak megfogalmazásánál az egyenlet még nem tartalmazza az elektromos tagot, ezt majd később a lokális alaknál illesztjük be az elektro-mágneses mező munkájából származtatva, a mágneses tér munkáját elhanyagolva. A termodinamika I. főtétele rugalmas alakváltozás esetére.
dk + de = dq + dw , (2.8) ahol dk – az elemi mozgási energia növekmény, de – az elemi belső energia növekmény, dq – az elemi elnyelt hő növekmény, dw – az elemi külső munka növekmény.
10
2. Piezo-termo-elasztikus alapegyenletek
Ezt felírhatjuk az idő szerinti deriváltakra is: &, k& + e& = q& + w (2.9) ahol k – a fajlagos mozgási energia, e – a fajlagos belső energia, q – a fajlagos elnyelt hő, w – a fajlagos külső munka. A változók fölötti pont idő szerinti deriváltat jelent. Egy hővezető rugalmas testre a (2.9) egyenlet a következő alakban írható (globális alak): d v2 ρ e + dV = − ∫ q ⋅ ndA + ∫ ρRdV + ∫ v ⋅ σ ⋅ ndA + ∫ F ⋅ vdV , dt V∫ 2 A V A V
(2.10) ahol v – a sebesség, n – a dA felület normálvektora, q – a hőáram, R – a belső hőforrás erősség. A (2.10)-es főtétel tehát azt fejezi ki, hogy a belső energia és a mozgási energia összegének időbeli változása megegyezik a felületen átadott hő, a test belsejében keletkező hő, a test felületén működő erők és a térfogati erők teljesítményének összegével. Ez az egyenlet az ún. globális alak, azaz az egész rendszerre vonatkozik. A későbbiekben azonban célszerűbb számunkra ennek az egyenletnek a lokális alakja. A továbbiakban néhány átalakítást végzünk: d dρ de ρedV = ∫ ρ + e dV ∫ dt V dt dt V
(2.11) Ha ρ állandó, akkor a (2.11) a következő alakban írható: de d ρedV = ∫ ρ dV . ∫ dt dt V V
(2.12)
11
2. Piezo-termo-elasztikus alapegyenletek
A (2.10) bal oldalának másik tagja a következőképpen alakítható át, szintén a ρ állandó figyelembevételével: d 1 2 &&dV . ρv dV = ∫ ρv ⋅ v&dV = ∫ v ⋅ ρu dt V∫ 2 V V (2.13) A (2.5) mozgásegyenlet alapján:
∫ v ⋅ ρu&&dV = ∫ v ⋅ (σ ⋅ ∇ + F )dV .
V
V
(2.14) A (2.13)-at és (2.14)-et a (2.10)-be helyettesítve, a Gauss-Osztogradszkij-tételt felhasználva, majd átrendezve: de ∫V ρ dt + v ⋅ σ ⋅ ∇ + v ⋅ F − ρR − F ⋅ v + ∇ ⋅ q − ∇ ⋅ (v ⋅ σ ) dV = 0 , (2.15) azaz de
∫ ρ dt + v ⋅ σ ⋅ ∇ + v ⋅ F − ρR − F ⋅ v + ∇ ⋅ q − v ⋅ σ ⋅ ∇ − σ : ∇ o v dV =
V
de = ∫ ρ − ρR + ∇ ⋅ q − σ : ∇ o v dV = 0. dt V (2.16) A (2.16)-os egyenletből pedig megkapjuk a termodinamika I. főtételének lokális alakját: ρ
de = ρR - ∇ ⋅ q + σ : ∇ o v . dt
(2.17) Ez a lokális alak azonban csak hővezető rugalmas test esetén érvényes. Ezen kívül figyelembe kell még venni az elektromos erők munkáját is. Az elektromágneses konstitutív egyenletek a következők: D = εe E , (2.18) B = µmH ,
(2.19) ahol εe – a permittivitás, µm – a permeabilitás. Az elektromágneses erők fajlagos munkája: 1 1 w em = ( E ⋅ D + H ⋅ B ) = (ε e E ⋅ E + µ m H ⋅ H ) . 2 2 (2.20)
12
2. Piezo-termo-elasztikus alapegyenletek
Ebből az elektromos rész (a mágneses részt elhanyagolva): 1 w e = εe E 2 . 2 (2.21)
Ennek megfelelően a teljesítmény: dw e dE dD = εe E ⋅ = E⋅ . dt dt dt
(2.22) Így a (2.17)-es lokális alak a következő módon egészül ki: ρ
de dD = ρR - ∇ ⋅ q + σ : ∇ o v + E ⋅ . dt dt
(2.23) Az I. főtétel mellett szükséges még a termodinamika II. főtétele, vagy más néven a Clausius-Duhem egyenlőtlenség: ds q ρT − ρR + ∇ ⋅ q − ⋅ ∇T ≥ 0 , dt T (2.24) ahol s – az entrópia. A ρR szorzat helyére a (2.23)-et helyettesítve a (2.24) egyenlőtlenség a következő módon alakul: ds de dD q ρT − ρ + σ : ∇ o v + E ⋅ − ⋅ ∇T ≥ 0 dt dt dt T (2.25) Ezek után bevezetjük a szabadenergia függvényt (Helmholtz-függvény). Ez a függvény a belső energia és az entrópia és hőmérséklet szorzatának különbsége: Ψ = e − Ts . (2.26) Feltételezzük, hogy ez a szabadenergia függvény az alakváltozási tenzornak, az elektromos eltolás vektornak, valamint a hőmérsékletnek is függvénye: Ψ = Ψ ( ε , D, T ) . (2.27) A (2.25)-ben a belső energia helyére behelyettesíthető a (2.26), azaz a belső energia kifejezhető a szabad energia függvény segítségével és ezt helyettesíthetjük be a termodinamika II. főtételébe:
13
2. Piezo-termo-elasztikus alapegyenletek
ρT
ds dΨ ds dT dD q − ρ + ρT + ρs + σ : ∇ o v + E ⋅ − ⋅ ∇T ≥ 0 , dt dt dt dt dt T (2.28)
azaz dT dD q dΨ - ρ − ⋅ ∇T ≥ 0 . + ρs + σ :∇ov + E ⋅ dt dt T dt (2.29) Ugyanígy eliminálhatjuk a belső energiát a (2.23) egyenletből is: ds dT dD dΨ +T +s ρ . = ρR - ∇ ⋅ q + σ : ∇ o v + E ⋅ dt dt dt dt (2.30) A (2.27) miatt a szabadenergia függvény idő szerinti deriváltja a következőképpen írható: & = ∂Ψ ⋅ ε& + ∂Ψ T& + ∂Ψ ⋅ D& . Ψ ∂D ∂ε ∂T (2.31) A (2.31)-et behelyettesítve a (2.30)-ba és a (2.29)-be a következőket kapjuk: ∂Ψ ∂Ψ & ∂Ψ − σ : ε& + ρ + s T + ρ − E ⋅ D& + ρ(Ts& − R ) + ∇ ⋅ q = 0 ρ ∂ε ∂T ∂D (2.32) q dΨ ∂Ψ & ∂Ψ − σ : ε& + ρ + s T + ρ − E ⋅ D& + ⋅ ∇T ≤ 0 ρ T dε ∂T ∂D (2.33) A fenti egyenletnek illetve egyenlőtlenségnek minden esetben teljesülni kell, ε& , T& és D& tetszőleges értékeire. Mivel a zárójeles tagok nem függnek e deriváltaktól, ezért ezeknek azonosan nullának kell lenniük. Így a következő összefüggések kaphatók: σ =ρ
∂Ψ , ∂ε (2.34)
s=-
∂Ψ , ∂T (2.35)
E =ρ
∂Ψ , ∂D (2.36)
14
2. Piezo-termo-elasztikus alapegyenletek
∇ ⋅ q + ρ(Ts& - R) = 0 , (2.37) q ⋅ ∇T ≤ 0 , T (2.38) továbbá: Θ = T − T0 , (2.39) ahol T0 - a kezdeti hőmérséklet [K], Θ - a hőmérséklet-különbség. A (2.37) egyenletet szokás energia mérleg egyenletnek nevezni. Feltételezzük, hogy a hőmérséklet-különbség a kezdeti hőmérséklethez képest kicsi. Θ << 1 . T0 (2.40) A hőrugalmasságtanban szintén szerepel ez a hőmérséklet-különbség. Ott a szakirodalomban a hőmérséklet-különbségre vonatkozó szokásos megszorítás 100-200 K. A hőrugalmasságtani egyenletek levezetése hasonló a jelen levezetéshez. BALLA [22] dolgozatában kifejti, hogy az egyenletek levezetéséhez ennél szigorúbb megszorítás szükséges, ennek ellenére azonban az így levezetett egyenletek a kísérleti eredményekkel jó egyezést mutatnak. Ezek után tekintsük a szabad energia függvény Taylor-sorba fejtését a (0,0,T0) pont körül a másodfokú tagokig! Ahhoz, hogy a sorba fejtés során a másodrendű tagoknál magasabb rendű tagokat elhanyagolhassuk, szintén szükséges, hogy a hőmérsékletkülönbség, az alakváltozási tenzor elemei és az elektromos eltolás vektor elemei kicsik legyenek. Ez a szakirodalomban megadottnál szigorúbb megkötést jelent a hőmérsékletkülönbségre. Sorba fejtve:
ρΨ (ε, D, T0 + Θ) = ρΨ (0,0, T0 ) + ρ
∂Ψ ∂Ψ ∂Ψ ⋅D+ρ :ε+ρ ∂ε 0 ∂D 0 ∂T
Θ+ T0
∂ 2Ψ ∂ 2Ψ 1 ∂ 2Ψ 2 2 + ρ 2 : ε + ρ 2 : D + ρ 2 Θ 2 + 2 ∂ε 0 ∂D 0 ∂T 0 2ρ
∂ 2Ψ ∂ 2Ψ ∂ 2Ψ ⋅ DΘ . : ε ⋅ D + 2ρ : εΘ + 2ρ ∂ε ⋅ ∂D 0 ∂ε∂T 0 ∂D∂T 0 (2.41)
15
2. Piezo-termo-elasztikus alapegyenletek
A (2.41) egyenlet ε, D és Θ szerinti deriválásával előállítható a feszültségtenzor, az elektromos térerősség vektor és az entrópia a (2.34), (2.35) és (2.36) szerint. Tekintsük a test terheletlen, úgynevezett „természetes” állapotát! Ekkor ε = D = 0 és T=T0. Ebben az állapotban Ψ (0,0, T0 ) = 0 és a (2.41) jobb oldalán szereplő lineáris tagok 0-vá válnak, vagyis: σ =ρ
∂Ψ ∂Ψ ∂ 2Ψ ∂ 2Ψ ∂ 2Ψ =ρ +ρ 2 :ε +ρ ⋅D+ρ Θ ∂ε ∂ε 0 ∂ε ⋅ ∂D 0 ∂ε∂T 0 ∂ε 0
(2.42) E =ρ
∂Ψ ∂Ψ ∂ 2Ψ ∂ 2Ψ ∂ 2Ψ =ρ +ρ 2 ⋅D+ρ Θ :ε+ρ ∂D ∂D 0 ∂ε ⋅ ∂D 0 ∂D∂T 0 ∂D 0
(2.43) − ρs = ρ
∂Ψ ∂Ψ ∂ Ψ ∂ Ψ ∂ Ψ =ρ +ρ 2 Θ+ρ ⋅D+ρ :ε ∂Θ ∂T 0 ∂D∂T 0 ∂ε∂T 0 ∂T 0 2
2
2
(2.44) A „természetes” állapotban ε = D = 0 , Θ = 0 valamint σ = 0 , E = 0 és s = 0 . Ebből következően az alábbi konstansoknak is rendre nullának kell lennie: ∂Ψ ∂Ψ ∂Ψ ρ =0, ρ =0, ρ = 0. ∂ε 0 ∂D 0 ∂T 0
A megmaradó együtthatókra vezessük be a következő jelöléseket: 1 ∂ 2Ψ ρ =c, 2 ∂ε 2 0 ρ
∂ 2Ψ ~ = −h , ∂ε∂D 0
1 ∂ 2Ψ ~, ρ 2 =m 2 ∂D 0 ρ
∂ 2Ψ = −β , ∂ε∂T 0
1 ∂ 2Ψ ρ = −a , 2 ∂T 2 0 ρ
∂ 2Ψ = − ~p . ∂D∂T 0
Ezekkel az együtthatókkal a konstitutív egyenletek a következő formában adódnak: ~ σ = c : ε − h ⋅ D − βΘ , ~ ~ ⋅ D − ~pΘ , E = −h : ε + m
− ρs = − β : ε − ~p ⋅ D − aΘ .
(2.45) (2.46)
(2.47) A szakirodalomban szokás az elektromos térerősség vektor helyett az elektromos eltolás vektort kifejezni. A végeselemes egyenletekhez is célszerűbb ez a forma. Így a konstitutív egyenletek a következő alakban írhatók:
16
2. Piezo-termo-elasztikus alapegyenletek
σ = c : ε − h ⋅ E − βΘ (2.48)
D = h : ε + m ⋅ E + pΘ (2.49) − ρs = − β : ε − p ⋅ E − aΘ (2.50) A mátrixok rendjének jobb áttekinthetősége miatt a konstitutív egyenletek indexes alakban a következők:
σ ij = c ijkl ε kl − h ijk E k − β ij Θ , (2.51)
D = h ε kl + m Ek + p Θ , i
ikl
ik
i
(2.52) − ρs = − β ε ij − p Ei − aΘ . ij
i
(2.53) Az együttható mátrixok szimmetria tulajdonságai a következők:
c ijkl = c jikl = c ijlk = c klij , (2.54)
m =m , ik
ki
(2.55)
β =β , ij
ji
(2.56)
h
ikl
=h . ilk
(2.57) Jól látható a (2.48) egyenletből, hogy amennyiben az elektromos térerősség vektor és a hőmérséklet-különbség nulla, akkor a Hooke-törvényhez jutunk. Ha csupán az elektromos térerősség vektort tekintjük nullának, akkor a Duhamel-Neumann anyagtörvényt kapjuk. Hasonlóan a (2.49) egyenlet esetén, ha az alakváltozási tenzor és a hőmérséklet-különbség nulla, akkor az elektromosságtanból jól ismert lineáris anyagtörvény áll elő. A fentiekhez teljesen hasonló módon juthatunk el a hőáramra vonatkozó egyenlethez, a Fourier-féle hővezetési törvényhez.
q = − k ⋅ ∇T . (2.58) A piezo-termo-elasztikus mező egyenleteket úgy kaphatjuk meg, hogy a kapott konstitutív egyenleteket behelyettesítjük a (2.1)-(2.4) Maxwell egyenletekbe, a (2.5) mozgásegyenletbe és a (2.37) energia mérleg egyenletbe. Feltételezhetjük azt, hogy a mágneses tér elhanyagolható, vagy legalábbis időben állandó. Ezzel a (2.2) tovább alakítható:
17
2. Piezo-termo-elasztikus alapegyenletek
∇× E = 0 .
(2.59) A (2.59) egyenlet abban az esetben teljesül, ha E = −∇Φ , ahol Φ elektrosztatikus potenciál. Ezek után az elektromos térerősség helyére mindig ezt helyettesíthetjük. A (2.49)-et a (2.4)-be helyettesítve adódik az első piezo-termo-elasztikus mezőegyenlet:
h ⋅ ∇ : ∇ o u − m ⋅ ∇ ⋅ ∇Φ + p ⋅ ∇Θ − γ = 0 . (2.60) A (2.5) mozgásegyenletbe helyettesítve a (2.48) összefüggést azt kapjuk, hogy
&& + F − β ⋅ ∇Θ + h ⋅ ∇ ⋅ ∇Φ = 0 . c ⋅ ∇ : ∇ o u − ρu (2.61) A (2.37)-be pedig a (2.50) összefüggésből az entrópiát és a (2.58)-as Fourier törvényt helyettesíthetjük be. Mivel Θ kicsi, ezért T helyére közelítésként T0-t írhatunk:
& ) − ρR ) = 0 . & + aΘ k ⋅ ∇ ⋅ ∇T - (T0 (β : ∇ o u& − p ⋅ ∇Φ (2.62) A (2.60)-(2.62) egyenletek tehát a piezo-termo-elasztikus mezőegyenletek, amelyek elsődleges változói az elmozdulás, a hőmérséklet-különbség és az elektromos potenciál. Konkrét feladatok esetén még további tagok is elhanyagolhatóak. Általában ilyen a belső hőforrás, vagy az elektromos töltéssűrűség.
2.2 Az általánosított lineáris piezo-termo-elasztikus egyenletek A Fourier-féle hővezetésből az következik, hogy a testben végtelen nagy a hőterjedési sebesség. Emiatt olyan általánosított elméleteket dolgoztak ki, amelyekkel a hővezetési egyenlet hiperbolikussá válik és így hullám-típusú hőáram leírására alkalmasak. Ilyen egyenlet kidolgozására több módszer is van. A szakirodalomban a leginkább elterjedt a Fourier-törvény CATTANEO [68] és VERNOTTE [69] által javasolt általánosítása. Az általánosított hiperbolikus hővezetési egyenlet a következő: ∂ 1 + τ q = −k ⋅ ∇T , ∂t (2.63) ahol τ – az ún. relaxációs idő, anyagjellemző. Ezzel a módosítással elérhető, hogy a hővezetési egyenlet már tartalmaz egy hőáramsebesség tagot is. A τ relaxációs idő meghatározásával, becslésével többen foglalkoztak (CHESTER [70], FRANCIS [71]), értéke 10-10s - 10-14s (gázok - fémek) között lehet. A relaxációs idő kis értéke miatt általában makroszkopikus méretekben nem
18
2. Piezo-termo-elasztikus alapegyenletek
feltétlen szükséges a hiperbolikus hővezetési egyenlet használata, azonban mikro- és nanométeres tartományban az esetek többségében már ez írja le megfelelően a hővezetési folyamatot [72], [73]. Az általános hővezetési egyenlet kapcsán a mai napig számos publikáció jelenik meg, újabb hőterjedési mechanizmusokat feltételezve a (2.63)-as egyenlethez képest még más tagokat építenek be. Ugyanakkor vannak, akik megkérdőjelezik a hiperbólikus hővezetési egyenlet létjogosultságát, vagy további problémákat vetnek fel (ZANCHINI [74]) vele kapcsolatban. Az általánosított Hővezetési egyenlettel kapcsolatos vitás kérdésekről SZEKERES [75] cikke ad áttekintést. Dolgozatomban az eredeti CATTANEO és VERNOTTE által javasolt (2.63)-as egyenletet használom. Azt, hogy mikronos mérettartományban már jelentősége lehet a hiperbólikus hővezetési egyenletnek, az a dimenziótlanított hővezetési egyenlet alapján látható be. Tekintsük a következő egydimenziós módosított hővezetési egyenletet: ∂ ∂T ~ ∂ 2 T + =a 2 , 1 τ ∂t ∂t ∂x
(2.64) ahol k ~ a – a hőmérséklet vezetési tényező, anyagjellemző ( ~ a= ). ρ⋅c A (2.64)-et tovább alakítva: ∂T ∂ 2T ∂ 2T +τ 2 =~ a 2 ∂t ∂t ∂x (2.65) Dimenziótlanítva a távolságot és a hőmérsékletet: ξ=
x , L (2.66)
ϑ=
T(x, t) − T∞ , T0 − T∞
(2.67) ahol L – a jellemző méret (pl. a fal vastagsága), T∞ - a végtelen távoli pont hőmérséklete, T0 - a fal hőmérséklete a kezdeti időpontban. Ezekkel a (2.65): ∂ 2ϑ 1 ∂ϑ 1 ∂ 2ϑ ( ) ( ) ⋅ (T0 − T∞ ) ⋅ − = ⋅ ⋅ − + ⋅ ⋅ T T T T τ 0 ∞ 0 ∞ ~ ~ a ∂t a ∂t 2 L2 ⋅ ∂ξ 2 (2.68)
19
2. Piezo-termo-elasztikus alapegyenletek
L2 ∂ϑ L2 ∂ 2 ϑ ∂ 2 ϑ ⋅ + τ⋅ ~ ⋅ 2 = 2 ~ a ∂t a ∂t ∂ξ (2.69) Bevezetve a Fourier-számot: F0 =
~ a ⋅t L2 (2.70)
A Fourier-számra fennáll: L2 dF0 = dt ~ a (2.71) Így a (2.69) a következő képpen írható: ~ a ∂ 2ϑ ∂ 2ϑ ∂ϑ + τ⋅ 2 ⋅ 2 = 2 ∂F0 ∂ξ L ∂F0 (2.72) Bevezetve a Π =
~ a ⋅τ jelölést: L ∂ 2 ϑ ∂ϑ Π + = ∇ 2ϑ 2 ∂F0 ∂F0 2
(2.73) Tekintsünk két egyszerű számpéldát! A hőmérséklet vezetési tényező általában 10 -10-5 m2/s nagyságrendű. A relaxációs idő nagyságrendje pedig a 10-10-10-14 s nagyságrendjébe esik. Azt vizsgálom, hogy a hiperbolikus tagnak milyen L jellemző méret esetén lesz 1%-4
os hatása ( Π 2 = 0.01 ). 1. példa: 2 −4 m ~ a ≈ 10 , τ = 10 −10 s , Π = 0.1 s
L = 10 −5 m 2. példa: 2 −5 m ~ a ≈ 10 , τ = 10 −10 s , Π = 0.1 s
L = 3,16 ⋅ 10 −6 m A jellemző méret a felvett anyagjellemzőkkel, tehát a mikronos mérettartományba esik. Vagyis megállapítható, hogy mikronos mérettartományban, illetve az alatt indokolt lehet a módosított hővezetési egyenlet használata, érdemes ezt az adott konkrét esetben mérlegelni.
20
2. Piezo-termo-elasztikus alapegyenletek
Hogyan befolyásolja ez a piezo-termo-elasztikus mezőegyenleteket? A Fourier-féle hővezetési egyenletet csak a (2.62)-es egyenletben használtuk fel, tehát az új hővezetési egyenlet is ezt az egyenletet befolyásolja. A (2.37)-be behelyettesítjük a (2.50) konstitutív egyenletet, a hőáram vektor esetén azonban nem használjuk a Fourier-törvényt. Így a (2.62)-as egyenlethez jutunk:
( (
)
)
& − ρR = 0 . & + aΘ − ∇ ⋅ q - T0 β : ∇ o u& − p ⋅ ∇Φ (2.74) A fenti egyenletet más alakban is írhatjuk az egyszerűbb áttekintés végett: & = 0. & − aT Θ − ∇ ⋅ q + ρR - T0 β : ∇ o u& + T0 p ⋅ ∇Φ 0 (2.75) Az aT0 szorzat helyére szokás bevezetni a ρc e szorzatot, ahol ce a fajhő (állandó alakváltozás mellett): & =0 & − ρc e Θ − ∇ ⋅ q + ρR - T0 β : ∇ o u& + T0 p ⋅ ∇Φ (2.76) ∂ Most alkalmazzuk az 1 + τ operátort mindkét oldalra: ∂t && + Θ & )= 0. && + ∇Φ & ) − ρc (τΘ && + ∇ o u& ) + T0 p ⋅ (τ∇Φ k ⋅ ∇ ⋅ ∇Θ + ρ(τR& + R ) - T0 β : (τ∇ o u e (2.77) Az általánosított hővezetési egyenlettel a piezo-termo-elasztikus mezőegyenletek végül is a (2.60), (2.61) és a (2.77) egyenletek.
2.3 Kezdeti és peremfeltételek A piezo-termo-elasztikus mezőegyenletek a fentiek alapján ismertek. Ahhoz, hogy a feladat teljes legyen, szükségesek még a kezdeti és peremfeltételek. A kezdeti és peremfeltételeket feloszthatjuk mechanikai, termikus és elektromos típusú mellékfeltételekre. Az alábbiakban közölt kezdeti és peremfeltételek mind a klasszikus, mind az általánosított hővezetési egyenlet esetén megfelelnek. Kezdeti feltételek: u( x,0) = u0 ( x ) ,
u& ( x,0) = u& 0 ( x ) ,
Θ( x ,0) = Θ 0 ( x ) ,
& ( x ,0) = Θ & ( x) , Θ 0
Φ( x,0) = Φ 0 ( x ) ,
& ( x,0) = Φ & 0 ( x) . Φ
21
2. Piezo-termo-elasztikus alapegyenletek
Peremfeltételek: A mechanikai peremfeltételeken belül megkülönböztetünk kinematikai és dinamikai peremfeltételeket. A kinematikai peremfeltétel azt jelenti, hogy az A felületű és V térfogatú test egy Au felületén adott az elmozdulás vektor, illetve egy Av felületén a sebességvektor. A dinamikai peremfeltétel pedig azt jeleneti, hogy egy Ap felületen ismert a külső erőrendszer felületre merőleges irányú sűrűség vektora. Mechanikai peremfeltételek: - kinematikai: ~ ( x, t ) , u=u -
x ∈ Au ,
~& ( x, t ) , u& = u
x ∈ Av .
dinamikai: σn = ~p( x, t ) ,
x ∈ Ap .
A termikus peremfeltételek esetében megkülönböztetünk első-, másod- és harmadfajú peremfeltételeket. Az elsőfajú peremfeltételnél a hőmérséklet különbséget ismerjük a test egy AT felületén. A másodfajú peremfeltétel a hőáramvektor normális komponensét írja elő a felület egy Aq részén. A harmadfajú peremfeltételnél pedig szintén a hőáram normális komponensét írjuk elő, de ebben az esetben ez konvektív hőátadás útján jön létre a test egy Ah felületén. Termikus peremfeltételek: - elsőfajú: ~ x ∈ AT . Θ = Θ( x , t ) , -
másodfajú: qn = ~ q ( x, t ) ,
x ∈ Aq .
harmadfajú: qn = h * (Θ − Θ ∞ ) ,
x ∈ Ah .
A harmadfajú peremfeltétel esetén h * – a hőátadási tényező, Θ ∞ – a környezet referencia hőmérsékletéből számított hőmérséklet különbség. Megjegyzendő, hogy a sugárzás útján átvitt hőmérséklet a harmadfajú peremfeltételhez hasonló módon adható meg az ún. sugárzási hőátadási tényező segítségével [22]. Az elektromos peremfeltételnél az AΦ felületen írhatjuk elő a potenciált. Elektromos peremfeltétel: ~ Φ = Φ( x, t ) ,
x ∈ AΦ
2. Piezo-termo-elasztikus alapegyenletek
22
Új tudományos eredmény: Levezettem a kapcsolt piezo-termo-elasztikus anyag- és mezőegyenleteket egy módosított (CATTANEO és VERNOTTE által javasolt), hiperbolikus hővezetési egyenlet felhasználásával. A Fourier-törvény alkalmazásával a kapcsolt piezo-termo-elasztikus egyenletek már ismertek. A hiperbolikus hővezetési egyenlet felhasználása azért fontos, mert a hőterjedéssel kapcsolatos jelenségek így nem végtelen nagy sebességgel terjednek, és ennek a mikronos mérettartományban nagy jelentősége van. Az egyenletek lineárisan rugalmas, homogén, anizotróp esetben érvényesek, kis alakváltozás mellett. A mezőegyenletek változói az elmozdulás, a hőmérséklet-különbség és az elektromos potenciál. Ezek lehetnek a végeselemes egyenletek változói is. A fejezet témájában megjelent publikációim: [P6], [P15], [P23]
3. Piezo-termo-elasztikus problémák végeselemes megfogalmazása
3.
PIEZO-TERMO-ELASZTIKUS
PROBLÉMÁK
23
VÉGESELEMES
MEGFOGALMAZÁSA
Az előző fejezetben levezetett egyenletek egzakt megoldása általános esetben nem állítható elő, sok speciális esetben rendkívül bonyolult. A szakirodalomban találhatunk néhány példát, ahol speciális geometriákat vizsgálnak (fél sík, végtelen henger, végtelen hosszú szalag), azonban a gyakorlatban előforduló szerkezetek egzakt vizsgálata reménytelennek látszik. Az ilyen feladatok esetében numerikus módszereket kell alkalmazni. Kapcsolt parciális differenciálegyenletek megoldására több numerikus módszer is kínálkozik. Napjainkban rendkívül elterjedt és hatékony módszer a végeselem módszer. Ebben a fejezetben a piezo-termo-elasztikus egyenletek végeselemes átírásával foglalkozom. Kapcsolt feladatok egyenleteit megoldó numerikus módszerekről BALLA [22] dolgozatában részletes áttekintést ad, így én csak a legszükségesebb tudnivalókat közlöm. A fizikai jelenségeket alapvetően kétféle módon írhatjuk le: vagy lokálisan vizsgálva a testet (mezőt), vagy globálisan. Az első esetben a jelenséget parciális differenciálegyenlet (rendszer) írja le, a második esetben integrál egyenlet (rendszer). Az integrál egyenlet variációs elvből származtatható, de a differenciálegyenletből is megalkotható. A variációs elvekből származtatott integrál egyenletek előnye, hogy az egyenletek mögött fizikai tartalom van, és az így kapott egyenletek megoldása a peremfeltételeket is kielégíti, azonban gondot okoz az, hogy számos jelenség leírására nem létezik ún. természetes variációs elv (pl. hővezetés). Ebben az esetben a differenciálegyenletből létrehozható ún. formális variációs elv, de ezek mögött nem feltétlenül létezik fizikai tartalom. A kétfajta leírásnak megfelelően a numerikus módszerek is két részre oszthatók. A differenciálegyenleteket véges differencia módszerekkel oldhatjuk meg, azaz a differenciálegyenletet diszkretizáljuk. Ezek alapvető hátránya, hogy csak az előre meghatározott pontokban és időpillanatokban kaphatunk közelítést, a pontok közötti eloszlásról nincs információ. Problémát jelenthetnek a bonyolult geometriájú szerkezetek, illetve a magas dimenzió szám is. Ezzel szemben az integrál egyenletek megoldására alkalmas módszerekkel (súlyozott maradványok módszere, végeselem módszer) nincs ilyen probléma. A végeselem módszer napjainkban rendkívül elterjedt és a súlyozott maradványok módszerével olyan feladatok megoldására is alkalmas, amelyekre nem léteznek variációs elvek [22].
3.1 Végeselemes egyenletek hely szerinti diszkretizációja A piezo-termo-elasztikus egyenletek végeselemes megfogalmazásához a Galjorkin módszer jól használható, mivel egyszerű és így az integrál egyenletek előállításához sem természetes, sem formális variációs elvekre nincs szükség.
3. Piezo-termo-elasztikus problémák végeselemes megfogalmazása
24
A végeselemes leíráshoz szükség van az egyenletek hely szerinti diszkretizációjára. Amennyiben dinamikai problémáról lenne szó, az idő szerinti diszkretizáció is szükséges. Az elmozdulás mezőt (u), a hőmérséklet különbséget (Θ) és az elektromos potenciált (Φ) tekintjük elsődleges változóknak. A formafüggvények segítségével a V térfogatú, A felületű test felosztása után ezek a változók közelítése a következő alakban írhatóak fel egy elemen belül: n
~ (r ) = N u (r , t )u e , u ∑ i i i =1
(3.1a) n
~ Θ(r ) = ∑ N iΘ (r , t )Θ ie , i =1
(3.1b) n
~ Φ(r ) = ∑ N iΦ (r , t )Φ ie , i =1
(3.1c) ahol ~ ~ ~ ,Θ u , Φ - a mezőfüggvények közelítései, N iu , N iΘ , N iΦ - formafüggvények, uie , Θ ie , Φ ie - a mezőfüggvények csomóponti értékei, n – egy elem csomópontszáma, t – idő, r - helyvektor. Szokás a hely szerinti és az idő szerinti diszkretizációt szétválasztani, ú.n. szemidiszkretizációt alkalmazni, mert ezáltal a végeselem hálózat dimenziója csökkenthető, és az időtartomány megváltozásakor nem szükséges új végeselem hálózatot generálni. Szemidiszkretizáció esetén azt feltételezzük, hogy a formafüggvények nem függnek az időtől, ugyanakkor a csomóponti változók időfüggők: n
~ (r ) = N u (r )ue ( t ) , u ∑ i i i =1
(3.2a) n
~ Θ(r ) = ∑ N iΘ (r )Θie ( t ) , i =1
(3.2b) n
~ Φ(r ) = ∑ N iΦ (r )Φ ie ( t ) . i =1
(3.2c)
3. Piezo-termo-elasztikus problémák végeselemes megfogalmazása
25
Az elsődleges változók idő szerinti deriváltjai is könnyen előállíthatóak. n
~& (r ) = N u (r )u& e ( t ) , u ∑ i i i =1
(3.3a) n ~& & e (t) , Θ(r ) = ∑ N iΘ (r )Θ i i =1
(3.3b) n
~& & e (t) , Φ(r ) = ∑ N iΦ (r )Φ i i =1
(3.3c) n
&~&(r ) = N u (r )u u ∑ i &&ie (t ) , i =1
(3.4a) n &~& && e ( t ) , Θ(r ) = ∑ N iΘ (r )Θ i i =1
(3.4b) n
&~& && ie ( t ) , Φ(r ) = ∑ N iΦ (r )Φ i =1
(3.4c) Az idő szerinti diszkretizációtól a következő alfejezetben lesz szó. Tekintsük hát először a hely szerinti diszkretizációt!
3.1.1 Hely szerinti diszkretizáció a klasszikus hővezetési egyenlet esetén
Rendezzük a csomóponti értékeket egy oszlopvektorba! (Az időfüggést a csomóponti értékeknél és a helyfüggést a formafüggvényeknél továbbiakban nem írom ki.) ue Y = Θ e , Φ e
(3.5) ahol u1e Θ1e Φ1e e e e u2 Θ2 Φ e e e ,Θ = , Φ = 2 . u = ... ... ... e e e u3 Θ 3 Φ 3
(3.6a,b,c)
26
3. Piezo-termo-elasztikus problémák végeselemes megfogalmazása
A formafüggvényeket a fenti sorrendnek megfelelően mátrixba rendezhetjük: N u N = 0 0
0 Θ
N 0
0 0 , N Φ
(3.7) ahol N1u Nu = 0 0
[
N = N Θ
[
Θ 1
N Φ = N1Φ
0 N1u
0 0
N u2 0
0 N u2
0 0
0
N1u
0
0
N u2
N
Θ 2
N Φ2
Θ n
... N un ... 0 ...
0
0 N un 0
0 0 , N un (3.8a)
]
... N , (3.8b)
]
... N Φn . (3.8c)
Tekintsük most a homogenizált mozgásegyenletet, hővezetési egyenletet és a Maxwellegyenleteket! Ezeknek az egyenleteknek a felhasználásával az előző fejezetben bemutatott – elsősorban a dinamikai, másod- és harmadfajú termikus – peremfeltételeket egyszerűen figyelembe lehet venni. A mágneses mező elhanyagolásával a (2.2) egyenlet a következő alakot ölti: ∇× E = 0 .
(3.9) A (3.9) egyenlet automatikusan teljesül, ha az elektromos térerősség vektor potenciálos: E = −∇Φ .
(3.10) A Galjorkin-módszer egy ú.n. súlyozott maradvány módszer. A súlyozott maradvány módszerek lényege, hogy a megoldandó, homogén alakra hozott differenciálegyenletbe az egyenletek megoldásának közelítését helyettesítjük, a közelítés miatt egy nem zérus maradványt (R) kapunk. Legyen egy parciális differenciálegyenlet P, melyet az x(r) függvény kielégít a V tartományon: P( x ) = 0 . x (r , a 1 , a 2 ,.., a n ) Legyen az x(r) függvény egy közelítése az ~
(3.11) függvény. Ekkor nyilván
P( ~ x ) = R (r, a1 , a 2 ,.., a n ) ≠ 0 ,
(3.12)
27
3. Piezo-termo-elasztikus problémák végeselemes megfogalmazása
ahol R (r, a1 , a 2 ,.., a n ) - a maradvány függvény. A súlyozott maradványok módszere szerint az a közelítés a legjobb, amelyiknél a maradvány súlyozott középértéke nulla.
∫ S(r )R (r , a , a 1
2
,.., a n )dV = 0 ,
V
(3.13) ahol S(r ) - a súlyfüggvény. Ahhoz, hogy az összes a i paramétert meghatározhassuk, n db egyenletre van szükség. Ezt az egyenletrendszert úgy állíthatjuk elő, hogy n db súlyfüggvényt választunk, melyeknek egymástól lineárisan függetleneknek kell lenni:
∫ S (r )R (r , a , a i
1
2
,.., a n )dV = 0 ,
i=1..n,
V
(3.14) ahol Si (r ) - egymástól lineárisan független súlyfüggvények. A Galjorkin-módszer esetében az Si(r) súlyfüggvények a formafüggvények:
∫ N (r )R (r , a , a i
1
2
,.., a n )dV = 0 ,
i=1..n,
V
(3.15) A Galjorkin-módszerhez tehát a diszkretizálandó egyenletekből fel kell építenünk egy maradványvektort minden végeselemre. A maradványvektor általában nem zérus, hiszen a változókat közelítjük, pontos értékét nem tudjuk, így a homogenizált egyenletek sem teljesülnek.
&~& − σ~ ⋅ ∇ − F Reu ρu ~& ~ Θ ~ & ~& − T p ⋅ ∇Φ Re = R e = q ⋅ ∇ + T0 aΘ - ρR + T0 β : ∇ o u 0 ≠ 0. ~ Φ R e D ⋅∇ − γ (3.16) ~ A fenti egyenletben σ~ , q~ és D az anyagtörvények segítségével a mezőfüggvényekből számított értékeket jelölik. A Galjorkin-módszert alkalmazva a következő integrál egyenlethez jutunk:
28
3. Piezo-termo-elasztikus problémák végeselemes megfogalmazása
N uT ⋅ R u W eu Θ T Θe Θ T W e = W e = ∫ N ⋅ Re dV = ∫ N R e dV . T Ve W eΦ V e N Φ R Φe
(3.17) A globális mátrix differenciálegyenletet az elemekre kapott súlyozott maradványokat összegezve kapjuk: m
W = ∑ We = 0 , e=1
(3.18) ahol m – a végeselemek száma. Mint azt a (3.16)-nél láthattuk, a piezo-termo-elasztikus egyenletek súlyozott maradványvektora három részből tevődik össze: a mozgásegyenletből, a hővezetési egyenletből és a Maxwell-egyenletekből származó részből. Amennyiben ezeket a súlyozott maradványvektorokat meghatározzuk, úgy a piezo-termo-elasztikus egyenleteket végeselemes számításra alkalmas formában kapjuk meg. A W eu komponens az elem csomópontszámával megegyező, egymáshoz hasonló módon előállítható alblokkból épül fel: Weu 1 u W u We = e 2 . ... u We n
(3.19) E vektornak az i-edik blokkja - egy 3x1-es vektor - a következő módon áll elő:
(
)
&~& − σ~ ⋅ ∇ − F dV . W eu i = ∫ N iu ρu ve
(3.20) A fenti alakot célszerű átalakítani és tagonként integrálni: ~ ~ ~ ∫ N (σ ⋅ ∇ )dV = ∫ (N σ )⋅ ∇dV − ∫ σ ⋅ (N ∇ ) dV . u i
Ve
u i
u i
Ve
e 14 4244 3 V1 44244 3
I
II
(3.21) A (3.21)-ben az I-jelű térfogati integrál tovább alakítható egy felületi intergrállá a Gauss-Osztrogradszkij tétel segítségével:
29
3. Piezo-termo-elasztikus problémák végeselemes megfogalmazása
∫ (N σ~ )⋅ ∇dV = ∫ (N σ~ )⋅ n dA = ∫ N (σ~ ⋅ n) dA . u i
u i
Ve
u i
Ae
Ae
(3.22) Ennek az átalakításnak az az előnye, hogy a felületen lévő peremfeltételeket könnyen figyelembe lehet venni (ld. előző fejezet 2.3-as pontja), mivel az n vektor a dA felület normálvektora: N u (σ~ ⋅ n ) dA = N u ~p dA ,
∫
∫
i
Ae
i
Ape
(3.23) ahol, Ape – az Ap felületnek a vizsgált elemre eső része (részfelülete). Nem veszem figyelembe a szomszédos elemek hatásait, mert az egyenletek összegzésénél, azaz az egész test vizsgálatánál ezek a hatások kiejtik egymást a hatás ellenhatás törvénye miatt. Ezen átalakítások után a Weu komponens a következő alakban írható fel: ~ &&dV − N u FdV − W eu i = ∫ σ~ ( N iu ∇)dV + ∫ N iu ρu ∫ i Ve
Ve
Ve
∫N
u i
~p dA .
A pe
(3.24) ~ Annyit kell még megtenni, hogy a (3.24)-be a σ helyére be kell helyettesítenünk a (2.48) anyagtörvényt, mivel a mezőváltozónk az elmozdulás és nem a feszültség. W eu i =
∫ (c : ε + h ⋅ ∇Φ − βΘ )( N ∇)dV + ∫ N ρu&&dV − ∫ N ~
~
~
u i
u i
Ve
Ve
~
Ve
u i
FdV −
∫N
u i
~p dA =
A pe
~ u ~ u u ~ u ~ u ∫ c : ε ( N i ∇)dV + ∫ h ⋅ ∇Φ( N i ∇)dV − ∫ βΘ( N i ∇)dV + ∫ N i ρu&&dV − ∫ N i FdV −
Ve
Ve
Ve
−
∫N
Ve
u i
Ve
~p dA .
A pe
(3.25) ~ Az ε már kifejezhető az elmozdulással a (2.7) alapján. Ezt az átalakítást azonban most nem hajtom végre, később az együttható mátrixok előállításánál még visszatérek rá. A W eΘ komponens szintén az elem csomópontszámával megegyező, egymáshoz hasonló módón előállítható alblokkból épül fel:
W eΘ
WeΘ1 Θ W = e 2 . ... Θ We n
(3.26)
30
3. Piezo-termo-elasztikus problémák végeselemes megfogalmazása
Ennek a vektornak az i-edik blokkja - egy skalár - a következő módon áll elő: ~& ~& ~& + T p ⋅ ∇Φ WeΘ i = ∫ N iΘ q~ ⋅ ∇ + T0 aΘ - ρR + T0 β : ∇ o u dV . 0 V e
(3.27) Ezt az integrált is tagonként vizsgáljuk: ~ ~ ∫ N (q ⋅ ∇ ) dV = ∫ (N q )⋅ ∇dV - ∫ N Θ i
Θ i
Ve
Θ i
∇ ⋅ q~dV .
Ve
e 14 4244 3 V1 42 4 43 4
I
II
(3.28) A Gauss-Osztrogradszkij tétel felhasználásával a (3.28) I-gyel jelölt tagja a következő módon tovább alakítható: ~ ~ ~ ∫ (N q )⋅ ∇dV = ∫ (N q )⋅ n dA = ∫ N (q ⋅ n) dA . Θ i
Θ i
Ve
Θ i
Ae
Ae
(3.29) A fenti alak az előző fejezet 2.3-as pontjában ismertetett másod- és harmadfajú peremfeltételekkel alakítható tovább: ~ ∫ N (q ⋅ n)dA = ∫ N Θ i
Ae
Θ i
(
)
~ h * Θ − Θ ∞ dA +
A he
∫N
~ q dA .
Θ i
A qe
(3.30) Ahe – az Ah felületnek, Aqe – az Aq felületnek a vizsgált elemre eső része (részfelülete). Csakúgy, mint a (3.23)-ban, itt sem veszem figyelembe a szomszédos elemek hatásait. A (2.60) Fourier törvény segítségével a (3.28) II-vel jelölt tagja is átalakítható:
∫N
Θ i
Ve
)( )
(
~ ∇ ⋅ q~dV = − ∫ k ⋅ N iΘ ∇ ⋅ Θ∇ dV , Ve
(3.31) így a (3.28) a következő alakban írható: ~ ∫ N (q ⋅ ∇ ) dV = ∫ N Θ i
Ve
A he
Θ i
(
)
~ h * Θ − Θ ∞ dA +
∫N
A qe
Θ i
(
)( )
~ ~ q dA + ∫ k ⋅ N iΘ ∇ ⋅ Θ∇ dV . Ve
(3.32) Behelyettesítés után a (3.27) az alábbi alakban írható fel:
31
3. Piezo-termo-elasztikus problémák végeselemes megfogalmazása
WeΘ i =
∫N
Θ i
(
)
~ h * Θ − Θ ∞ dA +
A he
∫N
Θ i
A qe
(
)( )
~ ~ q dA + ∫ k ⋅ N iΘ ∇ ⋅ Θ∇ dV + Ve
~& ~& ~& dV − N Θ T p ⋅ ∇Φ + ∫ N iΘ T0 aΘ dV + ∫ N iΘ T0 β : ∇ o u dV - ∫ N iΘ ρR dV . ∫ i 0 Ve
Ve
Ve
Ve
(3.33) A W
Φ e
komponens szintén az elem csomópontszámával megegyező, egymáshoz
hasonló módon előállítható alblokkból épül fel:
W eΦ
WeΦ1 Φ W = e 2 . ... Φ We n
(3.34) Ennek a vektornak az i-edik blokkja – amely a termikus részhez hasonlóan szintén egy skalár - a következő módon állítható elő:
(
)
~ WeΦ i = ∫ N iΦ D ⋅ ∇ − γ dV Ve
(3.35) A (2.49) segítségével ez az egyenlet átalakítható:
((
)
)
~ ~ ~ − m ⋅ ∇Φ + pΘ ⋅ ∇ − γ dV . WeΦ i = ∫ N iΦ h : ∇ o u Ve
(3.36) A (3.36)-et tagonként integrálva kapjuk a W
Φ e
i-edik blokkját:
~ ~ ~ ⋅ ∇dV − N Φ m ⋅ ∇Φ WeΦ i = ∫ N iΦ h : ∇ o u ⋅ ∇dV + ∫ N iΦ p∇ΘdV − ∫ N iΦ γdV . ∫ i Ve
Ve
Ve
Ve
(3.37) A (3.1), (3.2), (3.3) egyenleteket a (3.25), (3.33), (3.37) egyenletekbe helyettesítve a piezohőrugalmasságtan mátrix differenciálegyenlete megalkotható. A mezőváltozókat a formafüggvény mátrixszal (3.7) és a csomóponti értékekkel (3.5) kifejezhetjük. ~ Y = N ⋅Y e . (3.38) Minden integrált úgy fejezetem ki, hogy abban a mezőváltozók szerepeljenek, ezért a (3.38) behelyettesítése után a mezőváltozók csomóponti értékei kiemelhetőek.
&&e + K euu ⋅ u e + K euΘ ⋅ Θ e + K euΦ ⋅ Φ e − Feu , W eu = M euu ⋅ u (3.39)
32
3. Piezo-termo-elasztikus problémák végeselemes megfogalmazása ΘΘ & e ΘΦ & e ΘΘ e Θ &e WeΘ = LΘu e ⋅ u + L e ⋅ Θ + L e ⋅ Φ + K e ⋅ Θ − Fe ,
(3.40) e ΦΦ e Φ WeΦ = K eΦu ⋅ u e + K ΦΘ e ⋅ Θ + K e ⋅ Φ − Fe .
(3.41) A fenti egyenletekből az egyenletrendszer mátrixformába átírható: M euu 0 0
&&e 0 0 0 u && e Θu 0 0 ⋅ Θ + Le e & & 0 0 Φ 0
0 LΘΘ e 0
0 u& e K euu & e LΘΦ e ⋅ Θ + 0 & e K Φu 0 Φ e
K euΘ K ΘΘ e K ΦΘ e
K euΦ u e Feu 0 ⋅ Θ e − FeΘ = 0 . e FeΦ K ΦΦ e Φ
(3.42) Azaz a mátrix differenciálegyenlet az alábbi alakban is felírható: ~& ~ &~& M ⋅Y + L ⋅Y + K ⋅Y = F . (3.43)
3.1.2 Az együttható mátrixok elemei a klasszikus hővezetési egyenlet esetén Az együttható mátrixokban a csomóponti változók kiemelése után csak az anyagjellemzők és a formafüggvények illetve azok deriváltjai szerepelnek. Az együttható mátrixok elemeit célszerű úgy megadni, hogy számítógépes programozásra alkalmas formában szerepeljenek. A végeselem módszer esetében megszokott eljárás az anyagtörvényeket vektoros formában fogalmazzák meg a mátrixos helyett, kihasználva a feszültségtenzor szimmetriáját. σˆ = Cˆ ⋅ εˆ − eˆ T ⋅ E − βˆ Θ , (3.44)
D = eˆ ⋅ εˆ + η ⋅ E + pΘ , (3.45) − ρs = − βˆ T ⋅ εˆ − κ ⋅ E − λΘ ,
(3.46) ahol σˆ , εˆ - feszültség és alakváltozás vektor (6x1 méretű), E , D - villamos térerősség és elmozdulás vektor (3x1), Θ - hőmérséklet különbség, Cˆ , eˆ , βˆ , η , p - rugalmas (6x6), piezoelektromos (3x6), hőfeszültségi (6x1), dielektromos (3x3), piroelektromos (3x1) együttható mátrixok, ρ - sűrűség, s - entrópia, κ (1x3), λ - együtthatók.
3. Piezo-termo-elasztikus problémák végeselemes megfogalmazása
33
Az együttható mátrixok elemei a (3.25), (3.33) és (3.37) alapján adhatók meg. A merevségi mátrix elemei: A K euu blokk már jól ismert a rugalmasságtanból, és a végeselem módszerben gyakran előfordul:
K euu = ∫∫∫ B T ⋅ Cˆ ⋅ B dxdydz . (3.47) ~ A (3.25) egyenletnél emiatt nem részleteztem azt, hogy a c ⋅ ε hogyan alakítható tovább. A többi almátrix is hasonlóan állítható elő:
K euΘ = − ∫∫∫ B T ⋅ βˆ ⋅ N Θ dxdydz , (3.48) K euΦ = ∫∫∫ B T ⋅ eˆ T ⋅ B dxdydz , K
ΘΘ e
~ ~ = ∫∫∫ B T ⋅ k ⋅ B dxdydz + ∫∫∫ h k N Θ ⋅ N Θ dxdydz ,
(3.49) (3.50)
K
Φu e
= ∫∫∫ B ⋅ eˆ ⋅ B dxdydz , T
(3.51) K ΦΘ = ∫∫∫ B T ⋅ p ⋅ N Θ dxdydz , e (3.52) K
ΦΦ e
= − ∫∫∫ B ⋅ η ⋅ B dxdydz , T
(3.53) ahol ∂ ∂x 0 0 B=∂ ∂y 0 ∂ ∂z
0 ∂ ∂y 0 ∂ ∂x ∂ ∂z 0
0 0 ∂ ∂ ∂x ∂x ∂ ∂ ~ ∂ Θ ∂z N u , B N és B = N Φ . = ∂y ∂y 0 ∂ ∂ ∂ ∂z ∂z ∂y ∂ ∂x
34
3. Piezo-termo-elasztikus problémák végeselemes megfogalmazása
A csillapítási mátrix elemei: Θ ˆ T LΘu e = ∫∫∫ T0 β ⋅ B ⋅ N dxdydz ,
(3.54) Θ Θ LΘΘ e = ∫∫∫ T0 a N ⋅ N dxdydz ,
ΘΦ e
L
(
= − ∫∫∫ T0 B ⋅ p ⋅ B ⋅ N T
T
Θ
) dxdydz .
(3.55) (3.56)
A tömeg mátrix elemei: A tömeg mátrixnak mindössze egy eleme van: T
M euu = ∫∫∫ ρN u ⋅ N u dxdydz . (3.57) A terhelés vektor elemei: Feu = ∫∫∫ N u ⋅ F dxdydz + ∫ N u ⋅ pˆ dA , Ap
(3.58) F = ∫∫∫ ρRN dxdydz − ∫ N ~ q dA + ∫ h *Θ ∞ N Θ dA , Θ e
Θ
Θ
Aq
Ah
(3.59) F = ∫∫∫ γN dxdydz . Φ e
Φ
(3.60)
3.1.3 Hely szerinti diszkretizáció a módosított hővezetési egyenlet esetén A módosított hővezetési egyenlet esetén a W eu és W eΦ blokkok változatlanul maradnak, azonban a W eΘ i-edik blokkja a (3.27) helyett a következő módon állítható elő: ~& ∂ ~& ~& − T p ⋅ ∇Φ WeΘ i = ∫ N iΘ 1 + τ q~ ⋅ ∇ + T0 aΘ - ρR + T0 β : ∇ o u dV . 0 ∂t Ve (3.61)
35
3. Piezo-termo-elasztikus problémák végeselemes megfogalmazása
∂ A 1 + τ szorzó kiemelhető az integráljel elé, ezért a tagonkénti integrálás ugyanúgy ∂t végezhető el, mint a (3.27) esetében. Most azonban végre kell hajtani a beszorzást is. Figyelembe kell venni továbbá azt is, hogy most a Fourier törvény helyett a (2.65) módosított hővezetési egyenletet használjuk:
(
)
∂ ∂ ∂ ~ WeΘ i = 1 + τ ∫ N iΘ h * Θ − Θ ∞ dA + 1 + τ ∫ N iΘ ~ q dA − 1 + τ ∫ N iΘ q~ ⋅ ∇ dV + ∂t Ve ∂t Aqe ∂t A he
∂ ∂ ~& ~& dV + 1 + τ ∫ N iΘ T0 aΘ dV + 1 + τ ∫ N iΘ T0 β : ∇ o u ∂t Ve ∂t Ve
∂ ~& − 1 + τ ∫ N iΘ T0 p ⋅ ∇Φ dV ∂t Ve
∂ - 1 + τ ∫ N iΘ ρR dV . ∂t Ve (3.62) A beszorzás után a W
Θ e
i-edik blokkja, a (2.77) módosított hővezetési egyenlet
felhasználásával a következő módon írható: WeΘ i =
∫N
A he
(
Θ i
)
(
~ ~& ~& h * Θ + τΘ − Θ ∞ − τΘ ∞ dA +
)
(
~ ~ ~ ∫ N (q + τq& )dA + ∫ k ⋅ (N ∇ )⋅ (Θ∇ ) dV + Θ i
Θ i
A qe
Ve
)
(
)
~& ~ ~& ~ && && ~& + τ∇ o u ~ && dV − N Θ T p ⋅ ∇ Φ + ∫ N iΘ T0 a Θ + τΘ dV + ∫ N iΘ T0 β : ∇ o u + τΦ dV ∫ i 0 Ve
Ve
- ∫ N iΘ ρ(R + τR& ) dV .
Ve
Ve
(3.63) Ennek megfelelően a differenciálegyenlet-rendszer módosul. A (3.39) és (3.41) változatlan marad, míg a (3.40) helyett a következő egyenlethez jutunk: ΘΦ & e e Θ &e WeΘ = C eΘu ⋅ u& e + C ΘΘ ⋅ Φ + K ΘΘ e ⋅ Θ + Ce e ⋅ Θ − Fe .
(3.64) Ezzel a mátrix differenciálegyenlet a következő lesz: M euu Θu Me 0
0 M ΘΘ e 0
K euu + 0 K eΦu
&&e 0 0 u && e Θu M ΘΦ e ⋅ Θ + Le && e 0 0 Φ
K euΘ K ΘΘ e ΦΘ Ke
0 ( ΘΘ Le 0
0 u& e & e LΘΦ e ⋅ Θ + & e 0 Φ
K euΦ u e Feu ( 0 ⋅ Θ e − FeΘ = 0 . e FeΦ K ΦΦ e Φ (3.65)
36
3. Piezo-termo-elasztikus problémák végeselemes megfogalmazása
3.1.4 Az együttható mátrixok elemei a módosított hővezetési egyenlet esetén Csakúgy, mint a 3.1.2 pontban, az együttható mátrixokban a csomóponti változók kiemelése után csak az anyagjellemzők és a formafüggvények, illetve azok deriváltjai szerepelnek. Az együttható mátrix elemei hasonlóan adhatóak meg, mint a klasszikus hővezetési egyenlet esetén, ez esetben azonban a (3.25), (3.63) és (3.37) alapján. A merevségi mátrix elemei: A merevségi mátrix elemei nem változnak, ugyanazok lesznek, mint a 3.1.2 pontban. Ez természetes is, hiszen a hiperbolikus hővezetési egyenlet csak a rendszer dinamikai tulajdonságait változtatja meg. A módosított hővezetési egyenlettel tehát a tömegmátrix, a csillapítási mátrix és a terhelés vektor változik. A csillapítási mátrix elemei: A csillapítási mátrix elemei közül csak az LΘΘ e változik meg, a többi nem változik: ( Θ Θ Θ Θ * LΘΘ e = ∫∫∫ T0 aN ⋅ N dxdydz + ∫ τh N ⋅ N dA . Ah
(3.66) A tömeg mátrix elemei: A tömeg mátrixnál az M euu változatlan marad, viszont három újabb almátrix jelenik meg: T
M eΘu = ∫∫∫ τT0 βˆ ⋅ B ⋅ N Θ dxdydz ,
(3.67) M
ΘΘ e
= ∫∫∫ τT0 aN ⋅ N dxdydz , Θ
Θ
(
(3.68)
)
M ΘΦ = − ∫∫∫ τT0 B ⋅ pT ⋅ B T ⋅ N Θ dxdydz . e (3.69) A terhelés vektor elemei: Az elözőekhez hasonlóan csak FeΘ változik meg, a többi komponens nem változik:
(
)
( & )N Θ dA . FeΘ = ∫∫∫ ρ(R + τR& )N Θ dxdydz − ∫ N Θ ~ q + τ~ q& dA + ∫ h * (Θ ∞ + τΘ ∞ Aq
Ah
(3.70) A Galjorkin módszer nagy előnye, hogy bármilyen más kapcsolt egyenlet átírására is alkalmas és nem szükséges természetes variációs elvek megléte. Arra azonban érdemes
37
3. Piezo-termo-elasztikus problémák végeselemes megfogalmazása
odafigyelni, hogy olyan esetekben, ahol a feszültségtenzor nem szimmetrikus (mikropoláris testek) az anyagtörvényeknél szereplő mátrixok dimenziói megváltoznak.
3.2 Végeselemes egyenletek idő szerinti diszkretizációja A végeselemes egyenletek időbeli megoldására számos módszer áll rendelkezésünkre (Laplace-transzformáció, időbeli végeselemek, véges differencia sémák). Az irodalomban az egyik legelterjedtebben használt módszer a Newmark-féle véges differencia séma. A Newmark-sémának is számos alakja ismert, nagy előnye azonban az, hogy levezethető belőle egy egylépéses módszer. Ez azért előnyös, mert a számítás során könnyedén meg lehet változtatni az időlépés nagyságát, míg ez a kétlépéses sémák esetében nem igaz. A megoldandó egyenlet akár a klasszikus, akár a módosított hővezetési egyenlet esetén a következő. ~ ~& ~ && M ⋅Y + L ⋅Y + K ⋅Y = F . (3.43) ~ Az időt ∆t időlépésekre osztjuk, az i-edik lépésben a csomóponti változók értéke Yi . A Newmark-módszer általánosan a következő módon írható fel [22]: ~ ~& ~ && M ⋅ Yi + L ⋅ Yi + K ⋅ Yi = Fi , (3.71) ~ ~& ~ && M ⋅ Yi +1 + L ⋅ Yi +1 + K ⋅ Yi +1 = Fi +1 , ~ ~ ~& ∆t Yi +1 = Yi + ∆tYi + 2
2
(3.72)
((1 − 2η)Y~&& + 2ηY~&& ), i +1
i
(3.73)
)
(
~& ~& ~ ~ && && Yi +1 = Yi + ∆t (1 − ζ )Yi + ζYi +1 . (3.74) A (3.73)-ban és (3.74)-ben szereplő η a numerikus csillapítási (erősítési) paraméter, a ζ paraméter a gyorsulás változását szabályozza. Az η , ζ paraméter-párokra számos kombináció létezik. Az egyik leggyakrabban használt paraméterpár az η =1/4, ζ =1/2 (trapéz szabály). Bizonyítható, hogy ez adja a legpontosabb megoldást. ~ A fenti egyenletrendszert célszerű úgy átalakítani, hogy az Yi +1 kifejezhető legyen az iedik időlépésben kapott értékekkel [22]:
[M + ζ∆tL + η(∆t ) K ]⋅Y~ 2
i +1
[
]
~ ~& 2 = [M + ζ∆tL]⋅ Yi + ∆tM + (ζ − η)(∆t ) L ⋅ Yi +
1 1 && 2 3 ~ 2 + − η (∆t ) M + ζ − η (∆t ) L ⋅ Yi + η(∆t ) Fi +1 . 2 2 (3.75)
3. Piezo-termo-elasztikus problémák végeselemes megfogalmazása
38
~ ~ ~& && Amennyiben Yi +1 már ismert, a (3.73) és (3.74) egyenletek segítségével Yi +1 és Yi +1
már kiszámítható. A (3.75)-ös séma egylépéses, hiszen az i-edik lépésből számítható az i+1-edik. Az egylépéses módszer előnye, hogy az időlépés menet közben megváltoztatható. Erre azért van szükség, mert a gyors lefolyású tranziens jelenségekhez kisebb időlépés szüksége, míg az állandósult megoldáshoz már nagyobb időlépés is elegendő. A séma önindító is, csupán az első lépésben kell megadnunk a csomóponti változók és deriváltjaik értékeit. A fent ismertetett egylépéses módszer hátránya, hogy minden lépésben ki kell számítani a mezőváltozók deriváltjait is. Ha ezekre nincs szükségünk, akkor a feladat megoldásához nagyobb tároló területet használunk, mint amekkorára valójában szükség lenne. Azonban úgy gondolom, hogy a számítógépek mai teljesítménye és kapacitása mellett ez már nem jelent komoly problémát.
Új tudományos eredmény: Galjorkin módszerének felhasználásával előállítottam a kapcsolt piezo-termo-elasztikus egyenletek végeselemes diszkretizációját: a) a klasszikus, Fourier-féle hővezetési egyenlet alkalmazásával. b) a módosított, hiperbolikus hővezetési egyenlet felhasználásával. Mindkét esetre levezettem a diszkretizált mátrix differenciál egyenletet, amely alapján végeselemes algoritmus készíthető.
A piezo-termo-elasztikus egyenletek végeselemes diszkretizációjával előállított végeselemes mátrix-differenciálegyenletre vonatkozóan az irodalomban már található egykét példa, azonban ezek különbözőek. Az általam Galjorkin módszerével előállított egyenlet a [20] és [21] cikkek eredményeit támasztja alá. A hiperbolikus hővezetési egyenlet felhasználásával levezetett mátrixdifferenciálegyenlet, mások által még nem publikált. A fejezet témájában megjelent publikációim: [P6], [P15], [P23]
4. Végeselem program piezo-termo-elasztikus feladatokhoz (PiTE)
4. VÉGESELEM (PITE)
39
PROGRAM PIEZO-TERMO-ELASZTIKUS FELADATOKHOZ
A végeselemes egyenletek levezetése alapján egy végeselem programot készítettem. A végeselem program FORTRAN programnyelven (Fortran 90) készült és Smith és Griffiths [66] által publikált végeselem program moduljait is használja. Ezek a modulok általában a mátrix műveletek elvégzését teszik lehetővé (invertálás, determinánsszámítás), azonban a program lényegi része, az új elemtípus létrehozása és az együttható mátrixok összeállítása saját eredmény. Az előző fejezetben ismertetett végeselemes mátrix differenciálegyenletet statikus esetben, két lépésben (szekvenciálisan) is megoldható, hiszen a termikus rész külön leválasztható, a termikus kapcsoló tagok a piezoelektromos egyenlet terhelési oldalán vehetők figyelembe. Az így megmaradó piezoelektromos rész merevségi mátrixa szimmetrikus. Az általam írt programban, az egyenleteket egy lépésben oldom meg. Igaz ugyan, hogy így a merevségi mátrix nem lesz szimmetrikus, de a későbbiekben bonyolultabb anyagtörvények esetén is használható marad a program.
4.1 Elemtípus A végeselem modell elkészítéséhez szükséges egy elemtípus, ami alkalmas az adott feladat megoldására. Különböző típusú feladatok esetén (sík, tér, hengerszimmetrikus), különböző elemeket használhatunk, illetve ezen feladatokon belül is többféle elemtípus létezhet. Térbeli piezotermo-elasztikus feladatok megoldásához olyan elemre van szükség, amelynek öt szabadságfoka van csomópontonként, a három 4.1. ábra: 8 csomópontos piezo-termoelmozdulás komponens, a elasztikus elem hőmérséklet különbség és az elektromos potenciál. Az általam választott elemtípus a 8 csomópontos tégla elem (4.1. ábra). Ennek előnye az, hogy térbeli, tehát nem kell a feladatot leszűkíteni. A tégla elem választásának másik előnye, hogy példafeladatok könnyen felépíthetőek belőle. Az elem megadásához szükségesek a formafüggvények. Az elmozdulás mezőhöz használt formafüggvények ugyanazok lehetnek, mint az egyszerű 8 csomópontos tégla elem esetén. A hőmérséklet különbség és az elektromos potenciál esetében is a formafüggvényeknek ugyanazokat a kritériumokat kell teljesíteni, amelyeket az
4. Végeselem program piezo-termo-elasztikus feladatokhoz (PiTE)
40
elmozdulás mező esetén. Ezért a hőmérséklet különbség és az elektromos potenciál estében is ugyanazokat a formafüggvényeket használom, mint az elmozdulás mezőnél:
N iu = N iΘ = N iΦ , (4.1) 1 8 (1 − ξ )(1 − η)(1 − ζ ) 1 (1 − ξ )(1 − η)(1 + ζ ) 8 1 (1 + ξ )(1 − η)(1 + ζ ) 8 1 (1 + ξ )(1 − η)(1 − ζ ) . Nu = NΘ = NΦ = 8 1 (1 − ξ )(1 + η)(1 − ζ ) 8 1 (1 − ξ )(1 + η)(1 + ζ ) 8 1 (1 + ξ )(1 + η)(1 + ζ ) 8 1 8 (1 + ξ )(1 + η)(1 − ζ )
(4.2) A formafüggvények ilyen módon történő választása a formafüggvények deriváltjainak előállítását is megkönnyíti, azaz a hőmérséklet különbség és az elektromos potenciál esetében már nem szükséges újabb deriválást végezni, hiszen a formafüggvények deriváltjai az elmozdulás mezőnél már adottak. Az együttható mátrixok kompozíciójánál azonban figyelni kell, hogy a megfelelő deriváltak szerepeljenek a mátrix szorzatoknál.
4.2 Teszt példák Fortran nyelven egy végeselem programot készítettem, amely képes statikai piezotermo-elasztikus feladatok megoldására. Ellenőrizni kell, hogy a program megfelelően működik, ezért egy egyszerű teszt feladatott késztettem. A már ismert problémák ismeretében a programnak megfelelően kell működnie, vagyis a termikus és elektromos terhelések nélkül vissza kell adnia SMITH és GRIFFITHS eredményeit [66]. Az eredményeket ellenőriztem a COSMOS/M 2.0 és az ANSYS 8.0 programcsomagokkal is. Ezt a tesztet a program jól teljesítette (ld. 4.1-4.3. táblázatok). A második teszt a termomechanikai problémák megoldása. Ezeket az eredményeket csak a COSMOS/M 2.0 és ANSYS 8.0 programokkal tudtam ellenőrizni, a kapott értékek itt is jó egyezést mutatnak (ld. 4.4-4.7. táblázatok). Teljes piezo-termo-elasztikus probléma ellenőrzése lenne a legjobb teszt feladat, azonban ez meglehetősen bonyolult feladat. Kétféle módszer áll rendelkezésre: mérés vagy irodalmi példákkal történő összehasonlítás. A mérések
4. Végeselem program piezo-termo-elasztikus feladatokhoz (PiTE)
41
költségesek, hiszen a felhasznált anyag nem általános gépészeti anyag, számos anyagjellemző mérésre lenne szükség mielőtt a teszt feladat ellenőrzésére sor kerülhetne és a mikronos mérettartomány speciális, nagy pontosságú mérőberendezéseket igényelne. Az irodalomban kevés példa található, ezek általában speciálisak, és a megjelent cikkek alapján nem reprodukálhatóak.
4.2.1 Mechanikai analízis A modellt a 4.2. ábra szemlélteti. Az 1,4,5,8 csomópontok rögzítettek, tehát az 1,4,5,8 csomópontok által határolt lapnál van befogva. Terhelést a 11 és 12 csomópontoknál adjuk a testre. A terhelésnek erő, hő és villamos összetevői vannak. A test anyaga kadmiumszelenid, amelynek anyagjellemzői – irodalmi adatok alapján – a következők [6], [8]: C11 C 12 C Cˆ = 13 0 0 0 0 eˆ = 0 e 31 η11 η = 0 0
C12
C13
0
0
C 22
C 23
0
0
C 23
C 33
0
0
0 0
0 0
C 44 0
0 C 55
0
0
0
0
0
0
0
e15
0 e 32
0 e 33
e 24 0
0 0
0 η 22 0
βˆ = [β1 β 2
0 0 0 , 0 0 C 66
0 0 , 0
0 0 , η 33 β3
0 0 0] , T
p = [0 0 p 3 ] , T
C11=C22=74,1 GPa, C33=83,6 GPa, C12=45,2 GPa, C13=C23=39,3 GPa, C44=C55=13,17 GPa, C66=14,45 GPa, β1=β2=621000 N/mK, β3=551000 N/mK,
η1=η2=82,5 x 10-12 C2/mN, η3=90,2 x 10-12 C2/mN, p3=-2,94 x 10-6 C/m2K, e1=e2=-0,160 C/m2, e3=0,347 C/m2, e5=-0,138 C/m2, k=9 W/mK.
4. Végeselem program piezo-termo-elasztikus feladatokhoz (PiTE)
42
4.2. ábra: A teszt modell A terhelések: a 11 és 12 csomópontokon -1000 N nagyságú z irányú erő a mechanikai és a piezoelektromos esetben illetve a termikus esetben emellett ezekben a csomópontokban a hőmérséklet 30 oC. A referencia hőmérséklet 0 oC. A modell mérete: 5x5x10 mm3 Első esetben csak mechanikai terhelést adok a testre. Az általam írt programban a termikus és kapcsolt mezőkre vonatkozó anyagjellemzők befolyásolhatják a végeredményt, ezért ezeket a mátrixokat nullával kellene egyenlővé tenni. Azonban így a rendszer merevségi mátrixa szinguláris lenne, ezért egy-egy nem nulla elemet ( 10 −6 ) a mátrixokban hagyok. A Smith és Griffiths által közétett program csak izotróp anyagok modellezését teszi lehetővé. Ezért az ő programjukkal való ellenőrzésre a fenti Cˆ mátrix elemeiből „egyenértékű” izotróp anyagjellemzőt kell származtatni. Az így kapott izotróp anyagjellemzőket csak Smith és Griffiths programjában fogom használni a saját programomban a COSMOS-ban és az ANSYS-ban az anizotróp (ortotróp) jellemzőket adom meg. Az izotróp anyagjellemzők számításánál több lehetőségünk is van, hisz a rugalmassági mátrix bármely két elemét felhasználhatjuk egy ν és egy E érték meghatározásához. Ezek az értékek egymástól kisebb-nagyobb mértékben eltérőek lesznek. A számítások során a C11 és C 66 elemeket választottam. Az E és a ν értékét a következő egyenletekből kaphatjuk meg: C11 =
(1 − ν ) ⋅ E , (1 - 2ν )⋅ (1 + ν ) (4.3)
43
4. Végeselem program piezo-termo-elasztikus feladatokhoz (PiTE)
C 66 =
E . 2 ⋅ (1 + ν ) (4.4)
E két egyenletből kapott eredmények: E ≈ 39,8 GPa , ν = 0,38 A kapott csomóponti elmozdulás eredményeket az 4.1-4.3. táblázatok szemléltetik. Csomópont 1 2 3 4 5 6 7 8 9 10 11 12
Saját program Smith & Griffiths Cosmos/M 2.0 0,000E+00 0,000E+00 0,000E+00 4,704E-05 4,682E-05 4,560E-05 4,704E-05 4,682E-05 4,560E-05 0,000E+00 0,000E+00 0,000E+00 0,000E+00 0,000E+00 0,000E+00 -4,668E-05 -4,642E-05 -4,506E-05 -4,668E-05 -4,642E-05 -4,506E-05 0,000E+00 0,000E+00 0,000E+00 6,676E-05 6,564E-05 6,158E-05 6,676E-05 6,564E-05 6,158E-05 -6,897E-05 -6,939E-05 -6,598E-05 -6,897E-05 -6,939E-05 -6,598E-05 4.1. táblázat: x irányú elmozdulás [m]
Ansys 8.0 0,000E+00 4,704E-05 4,704E-05 0,000E+00 0,000E+00 -4,668E-05 -4,668E-05 0,000E+00 6,676E-05 6,676E-05 -6,897E-05 -6,897E-05
Csomópont 1 2 3 4 5 6 7 8 9 10 11 12
Saját program Smith & Griffiths Cosmos/M 2.0 0,000E+00 0,000E+00 0,000E+00 1,165E-05 1,130E-05 8,554E-06 -1,165E-05 -1,130E-05 -8,554E-06 0,000E+00 0,000E+00 0,000E+00 0,000E+00 0,000E+00 0,000E+00 -1,266E-05 -1,260E-05 -9,559E-06 1,266E-05 1,260E-05 9,559E-06 0,000E+00 0,000E+00 0,000E+00 4,440E-06 5,487E-06 3,370E-06 -4,440E-06 -5,487E-06 -3,370E-06 -4,554E-07 5,766E-07 6,145E-07 4,554E-07 -5,762E-07 -6,145E-07 4.2. táblázat: y irányú elmozdulás [m]
Ansys 8.0 0,000E+00 1,165E-05 -1,165E-05 0,000E+00 0,000E+00 -1,266E-05 1,266E-05 0,000E+00 4,440E-06 -4,440E-06 -4,554E-07 4,554E-07
44
4. Végeselem program piezo-termo-elasztikus feladatokhoz (PiTE)
Csomópont 1 2 3 4 5 6 7 8 9 10 11 12
Saját program Smith & Griffiths Cosmos/M 2.0 0,000E+00 0,000E+00 0,000E+00 7,547E-05 7,534E-05 7,658E-05 7,547E-05 7,534E-05 7,658E-05 0,000E+00 0,000E+00 0,000E+00 0,000E+00 0,000E+00 0,000E+00 7,362E-05 7,338E-05 7,482E-05 7,362E-05 7,338E-05 7,482E-05 0,000E+00 0,000E+00 0,000E+00 2,119E-04 2,095E-04 2,090E-04 2,119E-04 2,095E-04 2,090E-04 2,220E-04 2,230E-04 2,214E-04 2,220E-04 2,230E-04 2,214E-04 4.3. táblázat: z irányú elmozdulás [m]
Ansys 8.0 0,000E+00 7,547E-05 7,547E-05 0,000E+00 0,000E+00 7,362E-05 7,362E-05 0,000E+00 2,119E-04 2,119E-04 2,220E-04 2,220E-04
A jellemző elmozdulás z irányú. Ezeket az értékeket grafikonban is szemléltetem a 4.3. ábrán.
Z irányú elmozdulás 2,500E-04 2,000E-04 Saját program Smith & Griffiths Cosmos/M 2.0 Ansys 8.0
[m]
1,500E-04 1,000E-04 5,000E-05 0,000E+00 0
1
2
3
4
5
6
7
8
9 10 11 12
Csomópont
4.3. ábra: z irányú csomóponti elmozdulások [m] A 4.3 ábrán látható, hogy a kapott eredmények gyakorlatilag egybeesnek.
4. Végeselem program piezo-termo-elasztikus feladatokhoz (PiTE)
45
A Smith és Griffiths által írt program a feszültségeket csak az integrálási pontokban számítja ki, míg a végeselem programok a csomóponti feszültségeket adják meg. Az általam írt programban mind az integrálási pontokban, mind a csomópontokban kiszámítom a feszültségeket. Az x irányú feszültség értékek a 4.4-4.5. táblázatban és a 4.4. ábrán láthatóak. Elem
Integrálási pont Saját program Smith & Griffiths 1 -3,024E+08 -3,034E+08 2 2,993E+08 2,999E+08 3 -3,024E+08 -3,034E+08 4 2,993E+08 2,999E+08 1 5 -3,726E+08 -3,727E+08 6 -3,726E+08 -3,727E+08 7 3,757E+08 3,762E+08 8 3,757E+08 3,762E+08 1 -1,186E+08 -1,113E+08 2 1,469E+08 1,541E+08 3 -1,186E+08 -1,113E+08 4 1,469E+08 1,541E+08 2 5 -8,841E+07 -9,625E+07 6 -8,841E+07 -9,625E+07 7 6,010E+07 5,347E+07 8 6,010E+07 5,347E+07 4.4. táblázat: x irányú feszültség az integrálási pontokban[Pa]
Csomópont 1 2 3 4 5 6 7 8 9 10 11 12
Saját program Cosmos/M 2.0 Ansys 8.0 6,972E+08 6,750E+08 6,972E+08 2,695E+08 3,059E+08 2,695E+08 2,695E+08 3,059E+08 2,695E+08 6,972E+08 6,750E+08 6,972E+08 -6,918E+08 -6,678E+08 -6,918E+08 -2,967E+08 -3,545E+08 -2,967E+08 -2,967E+08 -3,545E+08 -2,967E+08 -6,918E+08 -6,678E+08 -6,918E+08 2,916E+08 2,964E+08 2,916E+08 2,916E+08 2,964E+08 2,916E+08 -2,425E+08 -2,071E+08 -2,426E+08 -2,425E+08 -2,071E+08 -2,426E+08 4.5. táblázat: x irányú csomóponti feszültség [Pa]
46
4. Végeselem program piezo-termo-elasztikus feladatokhoz (PiTE)
X irányú feszültség 8,000E+08 6,000E+08 4,000E+08 [Pa]
2,000E+08 0,000E+00 -2,000E+08 0 1 2 3 4 5 6 7 8 9 10 11 12
Saját program Cosmos/M 2.0 Ansys 8.0
-4,000E+08 -6,000E+08 -8,000E+08 Csomópont
4.4. ábra: x irányú csomóponti feszültségek [Pa] A többi feszültség komponens is hasonló egyezést mutat. A fentiek alapján megállapítható, hogy a program a mechanikai analízist jól végezte.
4.2.2 Termomechanikai analízis Második esetben a hőterhelést is vizsgáljuk. Ebben az esetben már csak a saját program a COSMOS/M és az ANSYS eredményeit tudjuk összehasonlítani. Az elmozdulásokat a 4.6-4.8. táblázatok és a 4.5. ábra szemléltetik. Csomópont 1 2 3 4 5 6 7 8 9 10 11 12
Saját program Cosmos/M 2.0 0,000E+00 0,000E+00 4,718E-05 4,573E-05 4,718E-05 4,573E-05 0,000E+00 0,000E+00 0,000E+00 0,000E+00 -4,655E-05 -4,494E-05 -4,655E-05 -4,494E-05 0,000E+00 0,000E+00 6,719E-05 6,199E-05 6,719E-05 6,199E-05 -6,841E-05 -6,541E-05 -6,841E-05 -6,541E-05 4.6. táblázat: x irányú elmozdulás [m]
Ansys 8.0 0,000E+00 4,703E-05 4,703E-05 0,000E+00 0,000E+00 -4,650E-05 -4,650E-05 0,000E+00 6,670E-05 6,670E-05 -6,826E-05 -6,826E-05
47
4. Végeselem program piezo-termo-elasztikus feladatokhoz (PiTE)
Csomópont 1 2 3 4 5 6 7 8 9 10 11 12
Saját program Cosmos/M 2.0 0,000E+00 0,000E+00 1,151E-05 8,413E-06 -1,151E-05 -8,413E-06 0,000E+00 0,000E+00 0,000E+00 0,000E+00 -1,277E-05 -9,669E-06 1,277E-05 9,669E-06 0,000E+00 0,000E+00 4,268E-06 3,211E-06 -4,268E-06 -3,211E-06 -8,173E-07 2,561E-07 8,173E-07 -2,561E-07 4.7. táblázat: y irányú elmozdulás [m]
Ansys 8.0 0,000E+00 1,167E-05 -1,167E-05 0,000E+00 0,000E+00 -1,284E-05 1,284E-05 0,000E+00 4,476E-06 -4,476E-06 -8,409E-07 8,409E-07
Csomópont 1 2 3 4 5 6 7 8 9 10 11 12
Saját program Cosmos/M 2.0 0,000E+00 0,000E+00 7,540E-05 7,646E-05 7,540E-05 7,646E-05 0,000E+00 0,000E+00 0,000E+00 0,000E+00 7,369E-05 7,496E-05 7,369E-05 7,496E-05 0,000E+00 0,000E+00 2,117E-04 2,086E-04 2,117E-04 2,086E-04 2,221E-04 2,216E-04 2,221E-04 2,216E-04 4.8. táblázat: z irányú elmozdulás [m]
Ansys 8.0 0,000E+00 7,529E-05 7,529E-05 0,000E+00 0,000E+00 7,360E-05 7,360E-05 0,000E+00 2,111E-04 2,111E-04 2,216E-04 2,216E-04
48
4. Végeselem program piezo-termo-elasztikus feladatokhoz (PiTE)
Z irányú elmozdulás 2,500E-04 2,000E-04 1,500E-04 [m]
Saját program Cosmos/M 2.0 Ansys 8.0
1,000E-04 5,000E-05 0,000E+00 0 1 2 3 4 5 6 7 8 9 10 11 12 Csomópont
4.5. ábra: z irányú csomóponti elmozdulások [m] Az eredményül kapott csomóponti hőmérsékleteket a 4.9. táblázat és a 4.6. ábra tartalmazza. Csomópont 1 2 3 4 5 6 7 8 9 10 11 12
Saját program Cosmos/M 2.0 0,00 0,00 12,19 12,19 12,19 12,19 0,00 0,00 0,00 0,00 10,78 10,78 10,78 10,78 0,00 0,00 15,94 15,94 15,94 15,94 30,00 30,00 30,00 30,00 4.9. táblázat: hőmérséklet [oC]
Ansys 8.0 0,00 13,55 13,55 0,00 0,00 15,00 15,00 0,00 18,39 18,39 30,00 30,00
49
4. Végeselem program piezo-termo-elasztikus feladatokhoz (PiTE)
Hőmérséklet 35,00 30,00
[m]
25,00 Saját program Cosmos/M 2.0 Ansys 8.0
20,00 15,00 10,00 5,00 0,00 0
1
2
3
4
5
6
7
8
9 10 11 12
Csomópont
4.6. ábra: csomóponti hőmérsékletek [oC] Az x irányú feszültség értékeket a 4.10. táblázat mutatja. Csomópont 1 2 3 4 5 6 7 8 9 10 11 12
Saját program Cosmos/M 2.0 Ansys 8.0 6,976E+08 6,778E+08 6,970E+08 2,714E+08 3,050E+08 2,700E+08 2,714E+08 3,050E+08 2,700E+08 6,976E+08 6,778E+08 6,970E+08 -6,915E+08 -6,660E+08 -6,891E+08 -2,945E+08 -3,538E+08 -2,974E+08 -2,945E+08 -3,538E+08 -2,974E+08 -6,915E+08 -6,660E+08 -6,891E+08 2,931E+08 2,967E+08 2,929E+08 2,931E+08 2,967E+08 2,929E+08 -2,355E+08 -2,108E+08 -2,461E+08 -2,355E+08 -2,108E+08 -2,461E+08 4.10. táblázat: x irányú csomóponti feszültség [Pa]
A termomechanikai analízisek alapján termomechanikai része is megfelelően működik.
megállapítható,
hogy
a
program
4. Végeselem program piezo-termo-elasztikus feladatokhoz (PiTE)
50
4.2.3 Piezoelektromos analízis Harmadik esetben a piezoelektromos viselkedést vizsgáljuk. Ebben az esetben már csak a saját program és az ANSYS eredményeit tudjuk összehasonlítani. Az elmozdulásokat a 4.11.-4.13. táblázatok és 4.7. ábra szemléltetik.
Csomópont Saját program Ansys 8.0 1 0,000E+00 0,000E+00 2 4,704E-05 4,705E-05 3 4,704E-05 4,705E-05 4 0,000E+00 0,000E+00 5 0,000E+00 0,000E+00 6 -4,669E-05 -4,667E-05 7 -4,668E-05 -4,667E-05 8 0,000E+00 0,000E+00 9 6,668E-05 6,684E-05 10 6,671E-05 6,681E-05 11 -6,907E-05 -6,888E-05 12 -6,903E-05 -6,891E-05 4.11. táblázat: x irányú elmozdulás [m]
Csomópont Saját program Ansys 8.0 1 0,000E+00 0,000E+00 2 1,165E-05 1,166E-05 3 -1,166E-05 -1,165E-05 4 0,000E+00 0,000E+00 5 0,000E+00 0,000E+00 6 -1,267E-05 -1,265E-05 7 1,266E-05 1,266E-05 8 0,000E+00 0,000E+00 9 4,503E-06 4,380E-06 10 -4,549E-06 -4,338E-06 11 -4,084E-07 -5,001E-07 12 3,617E-07 5,421E-07 4.12. táblázat: y irányú elmozdulás [m]
4. Végeselem program piezo-termo-elasztikus feladatokhoz (PiTE)
51
Csomópont Saját program Ansys 8.0 1 0,000E+00 0,000E+00 2 7,545E-05 7,548E-05 3 7,547E-05 7,547E-05 4 0,000E+00 0,000E+00 5 0,000E+00 0,000E+00 6 7,364E-05 7,360E-05 7 7,362E-05 7,362E-05 8 0,000E+00 0,000E+00 9 2,117E-04 2,120E-04 10 2,118E-04 2,120E-04 11 2,222E-04 2,219E-04 12 2,221E-04 2,219E-04 4.13. táblázat: z irányú elmozdulás [m]
Z irányú elmozdulás 2,500E-04 2,000E-04 1,500E-04 [m]
Saját program Ansys 8.0
1,000E-04 5,000E-05 0,000E+00 0 1 2 3 4 5 6 7 8 9 10 11 12 Csomópont
4.7. ábra: z irányú csomóponti elmozdulások [m]
Az eredményül kapott csomóponti elektromos potenciál értékeket a 4.14. táblázat és 4.8.-4.9. ábrák tartalmazzák.
4. Végeselem program piezo-termo-elasztikus feladatokhoz (PiTE)
52
Csomópont Saját program Ansys 8.0 1 0,0 0,0 2 -411,7 -411,6 3 1380,0 1380,0 4 0,0 0,0 5 0,0 0,0 6 -1418,0 -1417,4 7 449,1 449,0 8 0,0 0,0 9 -12750,0 -12745,0 10 -22960,0 -22956,0 11 22880,0 22874,0 12 12830,0 12828,0 4.14. táblázat: elektromos feszültség [V]
Elektromos potenciál 30000,0 20000,0
[V]
10000,0 0,0 0
1
2
3
4
5
6
7
8
9 10 11 12
-10000,0 -20000,0 -30000,0 Csomópont
4.8. ábra: csomóponti elektromos potenciál [V]
Saját program Ansys 8.0
4. Végeselem program piezo-termo-elasztikus feladatokhoz (PiTE)
53
4.9. ábra: csomóponti elektromos potenciál ANSYS 8.0-val [V] Az x irányú feszültség értékeket a 4.15. táblázat mutatja. Csomópont Saját program Ansys 8.0 1 6,971E+08 6,973E+08 2 2,691E+08 2,698E+08 3 2,692E+08 2,698E+08 4 6,972E+08 6,972E+08 5 -6,919E+08 -6,916E+08 6 -2,971E+08 -2,964E+08 7 -2,971E+08 -2,964E+08 8 -6,918E+08 -6,917E+08 9 2,900E+08 2,907E+08 10 2,898E+08 2,910E+08 11 -2,440E+08 -2,435E+08 12 -2,442E+08 -2,433E+08 4.15. táblázat: x irányú csomóponti feszültség [Pa] A piezoelektromos modell alapján is megállapítható, hogy a program piezoelektromos része megfelelően működik.
4. Végeselem program piezo-termo-elasztikus feladatokhoz (PiTE)
54
4.3 Mintapélda
4.3.1 Mindkét végén befogott tartó Számos mikromechanikai struktúrában találunk mindkét végén befogott tartóként modellezhető elemeket (gyorsulásérzékelők, gázérzékelők, tapintásérzékelők stb.). A modell a 4.10. ábrán látható. A megfogott csomópontok a rúd egyik végén (jelölje ezt a véget „A”) az 1,5,25,4,8,26,37,43,49, másik végén („B”) pedig a 21,23,35,22,24,36,42,48,54. Az „A” végen hőmérséklet változást nem engedünk meg, a „B” végen pedig az elektromos potenciál nulla. Négy esetet vizsgálok. Az első esetben a rúd terhelése tisztán mechanikai, a másodikban hőterhelés is fellép: a „B” keresztmetszetet állandó 200 oC-os hőmérsékleten tarjuk. A harmadik esetben nem lesz hőterhelés, de lesz elektromos terhelés: 25 kV az „A” keresztmetszetben. Az utolsó esetben pedig mind a három terhelést alkalmazzuk a modellre. A négy terhelési esetben kapott eredmények grafikonokon mutatom be. A grafikonokon a terhelési eseteket rendre a „Mechanikai”, „Termomechanikai”, „Piezoelektromos” és „Piezo-termo-elasztikus” jelzők hivatkozzák. A legjellemzőbb érték az x irányú feszültség, ezért az eredmények közül csak ezt mutatom be.
4.10. ábra: A mintapélda
55
4. Végeselem program piezo-termo-elasztikus feladatokhoz (PiTE)
4.11. ábra: A mintapélda modellje A 4.12.-4.14. ábrák az x-irányú feszültség vastagság menti eloszlását mutatják az „A” keresztmetszetben. A „B” keresztmetszetről hasonló ábrák készíthetőek, a számértékek eltérésével. Az ábrákon látható, hogy a két szélső vonalon az elektromos potenciál hatására a feszültségeloszlás lineáris jellege megmarad, azonban a meredeksége megváltozik, egyik oldalon nő, másik oldalon csökken. A középső vonalon az elektromos potenciálnak nincs hatása. x irányú feszültség 1,5E+08
Feszültség [Pa]
1,0E+08 5,0E+07 0,0E+00 -5,0E+07
0
10e-6
-1,0E+08
20e-6 Mechanikai Termomechanikai Piezoelektromos Piezo-termo-elasztikus
-1,5E+08 -2,0E+08
Vastagság [m]
4.12. ábra: Feszültség eloszlás a vastagság mentén (1,4,37 csomópontok)
56
4. Végeselem program piezo-termo-elasztikus feladatokhoz (PiTE)
x irányú feszültség 1,5E+08
Feszültség [Pa]
1,0E+08 5,0E+07 0,0E+00
0
10e-6
20e-6
-5,0E+07 -1,0E+08 Mechanikai Termomechanikai
-1,5E+08
Piezoelektromos Piezo-termo-elasztikus
-2,0E+08
Vastagság [m]
4.13. ábra: Feszültség eloszlás a vastagság mentén (5,8,43 csomópontok) x irányú feszültség 1,5E+08
Feszültség [Pa]
1,0E+08 5,0E+07 0,0E+00
0
10e-6
20e-6
-5,0E+07 -1,0E+08
Mechanikai Termomechanikai Piezoelektromos Piezo-termo-elasztikus
-1,5E+08 -2,0E+08
Vastagság [m]
4.14. ábra: Feszültség eloszlás a vastagság mentén (25,26,49 csomópontok) A 4.15.-4.17. ábrák a feszültségeloszlást szemléltetik a szélesség mentén az „A” keresztmetszetben.
57
4. Végeselem program piezo-termo-elasztikus feladatokhoz (PiTE)
x irányú feszültség 0,0E+00 -2,0E+07
0
20e-6
Mechanikai
-4,0E+07
Feszültség [Pa]
40e-6
Termomechanikai
-6,0E+07
Piezoelektromos -8,0E+07
Piezo-termo-elasztikus
-1,0E+08 -1,2E+08 -1,4E+08 -1,6E+08 -1,8E+08
Szélesség [m]
4.15. ábra: Feszültség eloszlás a szélesség mentén (1,5,25 csomópontok) Hasonlóan az előző három ábrához, ezekből is készíthető jellegre hasonló, számszerűleg némileg eltérő a „B” keresztmetszetben. Ezeken az ábrákon is megfigyelhető az, hogy a középső pontok a „Mechanikai”-„Piezoelektromos” és a „Termomechanikai”-„Piezotermo-elasztikus” esetekben egybe esnek. Ugyanakkor a szélső szálakon jelentős mértékű eltérés látható. x irányú feszültség 5,0E+06 0,0E+00
Feszültség [Pa]
-5,0E+06
0
20e-6
40e-6
-1,0E+07
Mechanikai
-1,5E+07
Termomechanikai Piezoelektromos
-2,0E+07
Piezo-termo-elasztikus
-2,5E+07 -3,0E+07 -3,5E+07 -4,0E+07 -4,5E+07
Szélesség [m]
4.16. ábra: Feszültség eloszlás a szélesség mentén (4,8,26 csomópontok)
58
4. Végeselem program piezo-termo-elasztikus feladatokhoz (PiTE)
x irányú feszültség 1,6E+08 1,4E+08
Feszültség [Pa]
1,2E+08 1,0E+08 8,0E+07 6,0E+07
Mechanikai
4,0E+07
Termomechanikai Piezoelektromos
2,0E+07
Piezo-termo-elasztikus
0,0E+00 0
20e-6
40e-6
Szélesség [m]
4.17. ábra: Feszültség eloszlás a szélesség mentén (37,43,49 csomópontok) A 4.18.-4.20. ábrákon a feszültség hossz menti eloszlása látható. Az előzőekben megállapítottak itt is érvényesek. x irányú feszültség Mechanikai Termomechanikai
1,0E+08
Piezoelektromos Piezo-termo-elasztikus
Feszültség [Pa]
5,0E+07 0,0E+00
0
20e-6
40e-6
60e-6
80e-6
100e-6
-5,0E+07 -1,0E+08 -1,5E+08 -2,0E+08
Hossz [m]
4.18. ábra: Feszültség eloszlás a hossz mentén (1,2,9,13,17,21 csomópontok)
59
4. Végeselem program piezo-termo-elasztikus feladatokhoz (PiTE)
x irányú feszültség Mechanikai Termomechanikai Piezoelektromos Piezo-termo-elasztikus
1,0E+08
Feszültség [Pa]
5,0E+07 0,0E+00
0
20e-6
40e-6
60e-6
80e-6
100e-6
-5,0E+07 -1,0E+08 -1,5E+08 -2,0E+08
Hossz [m]
4.19. ábra: Feszültség eloszlás a hossz mentén (5,6,11,15,19,23 csomópontok) x irányú feszültség Mechanikai
1,5E+08
Termomechanikai Piezoelektromos
1,0E+08
Feszültség [Pa]
Piezo-termo-elasztikus 5,0E+07 0,0E+00 -5,0E+07
0
20e-6
40e-6
60e-6
80e-6
100e-6
-1,0E+08 -1,5E+08 -2,0E+08 -2,5E+08
Hossz [m]
4.20. ábra: Feszültség eloszlás a hossz mentén (25,27,29,31,33,35 csomópontok) A rúd közepén tapasztalható némi aszimmetria. Ennek az az oka, hogy a modell valójában nem szimmetrikus, a tiszta mechanikai esetben sem, hisz a rúd „A” végén a hőmérsékletet, a „B” végén a feszültséget rögzítettem 0-ra.
4. Végeselem program piezo-termo-elasztikus feladatokhoz (PiTE)
60
Megállapítható ezen a mintapéldán, hogy az egyes terhelési esetek minőségileg és mennyiségileg is eltérő eredményt adnak.
Új tudományos eredmény: A 2. fejezetben levezetett egyenletek alapján végeselem algoritmust és programot készítettem. A programban a mátrix egyenletet - a korábbi publikációkban közöltekkel ellentétben, ahol a termikus részt leválasztva két lépésben (ún. szekvenciálisan) oldják meg az egyenleteket - egy lépésben (ún. direkt módon) oldom meg.
A direkt módszer előnye, hogy bonyolultabb anyagtörvények figyelembe vételével is használható (ahol a termikus rész nem választható le), hátránya, hogy a merevségi mátrix nem szimmetrikus. A fejezet témájában megjelent publikációim: [P6], [P15], [P23]
5. Mikromechanikai szerkezetek kapcsolt modellezése
61
5. MIKROMECHANIKAI SZERKEZETEK KAPCSOLT MODELLEZÉSE Ebben a fejezetben már konkrét mikromechanikai szerkezetek modellezésével foglalkozom, amelyek közül néhány meg is valósult. Ezeket a modelleket a COSMOS/M 2.0 programmal készítettem. A most bemutatandó szerkezetek anyaga piezoelektromos viselkedést nem mutat, a következő esetekben tehát termomechanikai analízist végeztem.
5.1. Egy dimenziós kapacitív gyorsulásérzékelő A gyorsulás érzékelésére számos elv létezik (5.1. ábra [76]). A mikromechanikai érzékelők körében leginkább elterjedt a piezorezisztív és a kapacitív típus. A kapacitív típusú gyorsulásérzékelő esetében egy rugalmasan felfüggesztett tömeg az érzékelő elem. A felfüggesztett tömeg a kondenzátor mozgó fegyverzete, vele szemben helyezkedik el egy rögzített fegyverzet. A tömeg gyorsulása arányos a rugó megnyúlásával, azaz a tömeg elmozdulásával, ez az elmozdulás pedig kapacitásváltozást okoz.
5.1. ábra: Gyorsulásérzékelő típusok sematikus rajza: a, piezoelektromos b, piezorezisztív c, rezgő elemes d, kapacitív e, optikai f, mágneses A felfüggesztett tömeget általában szilícium egykristályból alakítják, az általam bemutatott érzékelőt anizotróp marási technológiával. Ezt a réteget két üveg vagy szilícium réteg közé helyezik, ezek egyikén vagy mindkettőn alakítják a rögzített fegyverzetet. Két fajta kapacitív gyorsulásérzékelőt vizsgálok: egy oldalon felfüggesztett (lengőkaros) és két oldalon felfüggesztett (híd) típusút. A középső szilícium réteg rajzát az 5.2. ábrán láthatjuk.
62
5. Mikromechanikai szerkezetek kapcsolt modellezése
a, b, 5.2. ábra: egy oldalon (a) és két oldalon (b) felfüggesztett érzékelő
5.1.1. Mechanikai analízis A végeselem modell elkészítésénél a gyorsulásérzékelő mindhárom rétegét modelleztem. A modellben kihasználtam a szimmetriát, azaz a szerkezetnek csupán a felét készítettem el, a modellezéshez 8 csomópontos SOLID elemeket használtam. A felhasznált anyagjellemzők: Szilícium: E Si = 130 GPa
Üveg (#7740 pyrex): E Üveg = 63 GPa
ν Si = 0,3
ν Üveg = 0,2 3
ρSi = 2300 kg/m
ρ Üveg = 2200 kg/m3
E a Young-modulusz, ν a Poisson-tényező, ρ a sűrűség.
Az egy oldalon felfüggesztett és a két oldalon felfüggesztett szerkezetnél is a két üveg réteg megegyezik, csupán a középső szilícium réteg különbözik. A szerkezetre adott gyorsulás 10 g volt (y irányban).
5. Mikromechanikai szerkezetek kapcsolt modellezése
5.3. ábra: A háromrétegű gyorsulásérzékelő végeselem modellje
5.4. ábra: A középső, szilícium réteg az egy oldalon felfüggesztett érzékelőnél
5.5. ábra: A középső, szilícium réteg a két oldalon felfüggesztett érzékelőnél
63
64
5. Mikromechanikai szerkezetek kapcsolt modellezése
5.1.2. Termomechanikai analízis A termomechanikai modellnél már hőterhelést is figyelembe vettem. A modell, a mechanikai terhelés ugyanazok, mint az előző pontban, azonban további anyagjellemzőkre van szükség. Szilícium: α Si = 2,6⋅10-6 1/K
Üveg (#7740 pyrex): α Üveg = 3,25⋅10-6 1/K
k Si = 150 W/mK
k Üveg = 1,12 W/mK
c Si = 0,7 kJ/kgK
c Üveg = 0,83 kJ/kgK
α a hőtágulási együttható, k a hővezetési tényező, c a fajhő. A mechanikai terhelés mellett a következő termikus peremfeltételeket alkalmaztam a modellre: • A modell alsó lapján konvekció van, a közeg hőmérséklete 100 oC, a hőátadási tényező 10 W/m2K; • A felső és oldallapokon szintén konvekcióval történik hőátadás, a környezet 20 o C-os, a hőátadási tényező 5 W/m2K.
5.1.3. A mechanikai és termomechanikai eredmények összehasonlítása A modell jellemző geometriai méreteit különböző intervallumokban változtattam. • Felfüggesztő híd hossza: 200 µm < L1 < 1200 µm, L2=3000 µm • Felfüggesztő híd szélessége: 200 µm < s1 < 1200 µm, s2=3000 µm • Felfüggesztő híd vastagsága: 10 µm < v1 < 110 µm, v2=400 µm A termikus peremfeltételek következtében, a szilícium jó hővezetési képességei miatt, a szilícium rétegben közel homogén 47 oC-os hőmérséklet alakult ki. Kapacitív gyorsulásérzékelők esetén mind a tömeg elmozdulása, mind a felfüggesztő hidakban kialakuló feszültségek ismerete szükséges. Előbbi a szerkezet érzékenységét befolyásolja, utóbbi pedig a szerkezet terhelhetősége miatt fontos. Az 5.6.-5.11. ábrák a felfüggesztett tömeg lehajlását mutatják hossz- (L2/L1) szélesség(s2/s1) és vastagságarány (v2/v1) függvényében. Az 5.6. és 5.7. ábrákon látható, hogy a hídhossz növelésével (L2/L1 csökken), a lehajlás mértéke is nő. A termomechanikai modellezés során kapott értékek nagyobbak, mint termikus terhelés nélkül. A két oldalon felfüggesztett szerkezet esetében ez az eltérés jelentősebb, mint az egy oldalon felfüggesztett érzékelőnél.
65
5. Mikromechanikai szerkezetek kapcsolt modellezése
Függőleges elmozdulás 1
3
5
7
9
11
13
15
0 -0,5
Elmozdulás [mikron]
-1 -1,5 -2 -2,5 -3 -3,5 -4
Mechanikai Termomechanikai
-4,5 -5
L2/L1
5.6. ábra: Függőleges elmozdulás a hídhossz függvényében, egy oldalon felfüggesztett érzékelő esetében. Függőleges elmozdulás 1
3
5
7
9
11
13
15
1,00E-02
Elmozdulás [mikron]
5,00E-03 0,00E+00 -5,00E-03 -1,00E-02 -1,50E-02 -2,00E-02 -2,50E-02 -3,00E-02 -3,50E-02
Mechanikai Termomechanikai
-4,00E-02 -4,50E-02
L2/L1
5.7. ábra: Függőleges elmozdulás a hídhossz függvényében, két oldalon felfüggesztett érzékelő esetében.
66
5. Mikromechanikai szerkezetek kapcsolt modellezése
Az 5.8. és 5.9. ábrákon azt mutatják, hogy a hídszélesség csökkentésével (s2/s1 nő), a lehajlás mértéke nő. A termikus terhelés a két oldalon felfüggesztett szerkezet esetén jelentős, míg az egy oldalon felfüggesztet érzékelőnél csekély hatású. Függőleges elmozdulás 1
3
5
7
9
11
13
15
0
Elmozdulás [mikron]
-1
Mechanikai Termomechanikai
-2 -3 -4 -5 -6 -7 -8 -9 -10
s2/s1
5.8. ábra: Függőleges elmozdulás a hídszélesség függvényében, egy oldalon felfüggesztett érzékelő esetében. Függőleges elmozdulás 1
3
5
7
9
11
13
15
2,00E-02
Elmozdulás [mikron]
1,00E-02 0,00E+00 -1,00E-02 -2,00E-02 -3,00E-02 -4,00E-02 -5,00E-02 -6,00E-02 -7,00E-02
Mechanikai Termomechanikai
-8,00E-02
s2/s1
5.9. ábra: Függőleges elmozdulás a hídszélesség függvényében, két oldalon felfüggesztett érzékelő esetében.
67
5. Mikromechanikai szerkezetek kapcsolt modellezése
Az 5.10. és 5.11. ábrákon azt szemléltetik, hogy a hídvastagság csökkentésével (v2/v1 nő), a lehajlás mértéke ismét nő. A termikus terhelés a két oldalon felfüggesztett szerkezet esetén nagyobb hatású, míg az egy oldalon felfüggesztet érzékelőnél szinte elhanyagolható. Függőleges elmozdulás 2
7
12
17
22
27
32
37
42
0
Elmozdulás [mikron]
-50 -100 -150 -200 -250 -300 -350
Mechanikai Termomechanikai
-400
v2/v1
5.10. ábra: Függőleges elmozdulás a hídvastagság függvényében, egy oldalon felfüggesztett érzékelő esetében. Függőleges elmozdulás 2
7
12
17
22
27
32
37
42
0,00E+00
Elmozdulás [mikron]
-2,00E-01 -4,00E-01 -6,00E-01 -8,00E-01 -1,00E+00 -1,20E+00
Mechanikai Termomechanikai
-1,40E+00
v2/v1
5.11. ábra: Függőleges elmozdulás a hídvastagság függvényében, két oldalon felfüggesztett érzékelő esetében.
68
5. Mikromechanikai szerkezetek kapcsolt modellezése
Az 5.12.-5.17. ábrák a feszültségek alakulását mutatják a felfüggesztő hidakban. Egyenértékű feszültség (HMH) 6
Feszültség [MPa]
5
4
3
2
1
Mechanikai Termomechanikai
0 1
3
5
7
9
11
13
15
L2/L1
5.12. ábra: HMH-egyenértékű feszültség a hídhossz függvényében, egy oldalon felfüggesztett érzékelő esetében. Egyenértékű feszültség (HMH) 25
Feszültség [MPa]
20
15
10
5
Mechanikai Termomechanikai
0 1
3
5
7
9
11
13
15
L2/L1
5.13. ábra: HMH-egyenértékű feszültség a hídhossz függvényében, két oldalon felfüggesztett érzékelő esetében.
69
5. Mikromechanikai szerkezetek kapcsolt modellezése
Egyenértékű feszültség (HMH) 16 14
Feszültség [MPa]
12 10 8 6 4
Mechanikai Termomechanikai
2 0 1
3
5
7
9
11
13
15
s2/s1
5.14. ábra: HMH-egyenértékű feszültség a hídszélesség függvényében, egy oldalon felfüggesztett érzékelő esetében. Egyenértékű feszültség (HMH) 9 8
Feszültség [MPa]
7 6 5 4 3 2
Mechanikai Termomechanikai
1 0 1
3
5
7
9
11
13
15
s2/s1
5.15. ábra: HMH-egyenértékű feszültség a hídszélesség függvényében, két oldalon felfüggesztett érzékelő esetében.
70
5. Mikromechanikai szerkezetek kapcsolt modellezése
Egyenértékű feszültség (HMH) 140 120
Feszültség [MPa]
100 80 60 40
Mechanikai Termomechanikai
20 0 2
7
12
17
22
27
32
37
42
v2/v1
5.16. ábra: HMH-egyenértékű feszültség a hídvastagság függvényében, egy oldalon felfüggesztett érzékelő esetében. Egyenértékű feszültség (HMH) 1,00E+01 9,00E+00
Feszültség [MPa]
8,00E+00 7,00E+00 6,00E+00 5,00E+00 4,00E+00 3,00E+00 2,00E+00
Mechanikai Termomechanikai
1,00E+00 0,00E+00 2
7
12
17
22
27
32
37
42
v2/v1
5.17. ábra: HMH-egyenértékű feszültség a hídvastagság függvényében, két oldalon felfüggesztett érzékelő esetében.
5. Mikromechanikai szerkezetek kapcsolt modellezése
71
Alapvetően mindegyik ábráról leolvashatjuk, hogy a hőterhelés hatására az egyenértékű feszültség megnő. Lényeges különbség van azonban az egy, illetve a két oldalon felfüggesztett szerkezet között. Míg az egy oldalon felfüggesztett érzékelő esetében ez a növekedés kismértékű, addig a két oldalon felfüggesztett esetében drasztikus változásokat láthatunk. Azt is meg kell jegyezni, hogy a legnagyobb egyenértékű feszültség az egy hidas érzékelőnél nem a felfüggesztő hídban, hanem a keretben ébredt. A geometriai peremfeltételek tehát nagyban befolyásolják a hőterhelés hatására kialakuló feszültségeket, és számos esetben ezek a feszültségek hatása jelentős a szerkezet működésére.
5.18. ábra: Hőfeszültségek a keretben [MPa]
5.2. Gázérzékelő A gázérzékelő eszközök szintén vékony hidakon felfüggesztett struktúrák. A felfüggesztő hidak szerepe azonban itt teljesen más, mint a gyorsulásérzékelők esetében. A gázérzékelők magas hőmérsékleten működő eszközök, a hatékony működéshez feltétlenül szükséges a jó hőszigetelés. Ezért az érzékelő rész a levegőben „függeszkedik”, hiszen a levegő hővezetési tulajdonsága rossz. A szerkezet felfűtése elektromosan történik, a hidak szerepe az elektromos kapcsolat biztosítása is. Az elektromos hozzávezetés platina vagy többkristályos szilíciummal történhet, a hőszigetelés biztosítása érdekében ezeket szilícium-nitriddel (SixNy) vonják be. A szerkezet terhelése kizárólag hőterhelés, de a geometriai peremfeltételek miatt a keletkező hőfeszültségek kulcsfontosságúak lehetnek. A felfüggesztett lapkára még egy katalizátor réteg is felkerül. A szerkezet úgy működik, hogy a felfűtött lapka felett áramló gáz a katalizátor jelenléte miatt kémiai reakcióban vesz
5. Mikromechanikai szerkezetek kapcsolt modellezése
72
részt (elég), ennek következtében a felfűtött lap hőmérséklete tovább emelkedik. Ezt a hőmérséklet változást pedig elektromos úton lehet érzékelni.
5.19. ábra: Négy- és kéthidas gázérzékelő szerkezet A következőkben a fenti két szerkezetet vizsgálom, más-más szempontból. A négyhidas szerkezetnél a peremfeltételek hatása, a kéthidasnál az elektromos hozzávezetés anyaga lesz a vizsgálatok fókuszában.
5.2.1. Négyhidas gázérzékelő szerkezet A gázérzékelő modellezésénél az alapvető probléma az, hogy a peremfeltételek pontosan nem ismertek. Ezért különböző modelleket készítettem, amelyek segítségével a peremfeltételek hatása figyelemmel kísérhető. A modellezés során a hőterjedést kizárólag a hővezetésből számítottam a hőátadást és a sugárzást elhanyagoltam. Ez általában megtehető az eszköz méretei miatt. Két különböző típusú peremfeltételt alkalmaztam a modellre: PF1: A fűtőlap alatti lemez és a felfüggesztő hidak külső végeinek hőmérsékletét szobahőmérsékletűnek (20oC) tételezem fel. PF2: A fűtőlap alatti lemez, a felfüggesztő hidak külső végeinek és a környező levegő hőmérsékletét is szobahőmérsékletűnek (20oC) tételezem fel.
A modellezés során a fűtőlap alatti légrés méretét is változtattam, hogy ennek hatása is figyelembe vehető legyen. A szimmetria miatt elegendő csupán a szerkezet negyedét modellezni. A modellben a levegőt is modelleztem, a levegő hővezetésének figyelembe vételéhez.
73
5. Mikromechanikai szerkezetek kapcsolt modellezése
5.20. ábra: A fűtő elem modellje
5.21. ábra: A teljes modell hálózás nélkül
5.22. ábra: A peremfeltételek szemléltetése (PF1 és PF2) A modellezéshez szükségesek a felhasznált anyagok anyagjellemzői. Mikromechanikai anyagok esetén általános probléma az anyagjellemzők kérdése. Az irodalomban megtalálható anyagjellemzők általában nagy mértékben eltérnek egymástól, saját mérés végrehajtása pedig költséges és a méretek miatt körültekintő előkészítést kíván. Az általam vizsgált esetben a következő anyagjellemzőket használtam [77]: E [GPa] K [W/mK] C [J/kgK] ρ [kg/m3] α [1/K] ν SixNy Pt Levegő
270 170 -
0,27 0,3
20 75 0,024
2,11⋅10-6 -6
8,9⋅10
710,6
3100
130 1000
21440 1,1
5.1. táblázat: A felhasznált anyagok anyagjellemzői: E - Young-modulusz, ν - Poisson tényező, K - hővezetési tényező, α - hőtágulási együttható, C - fajhő, ρ - sűrűség.
5. Mikromechanikai szerkezetek kapcsolt modellezése
74
A fűtőelem terhelése 10 mW körüli elektromos teljesítmény, ez, illetve ennek egy része alakul hővé. A modellezés során azt feltételeztem, hogy a teljes 10 mW teljesítmény hővé alakul. Ezt a teljesítményt a fűtött rész belső, 80x80 µm2-es tartományára alkalmaztam (1560 mW/mm2). A modell jellemző méretei: a felfüggesztő hidak hossza 100 µm, szélessége 20 µm, vastagsága 1,6 µm. A fűtőlap mérete 100x100 µm2.
5.23. ábra: A kialakuló hőmérséklet eloszlás az 1. peremfeltételek mellett [oC], 20µm-es légrés mellett
5.24. ábra: A kialakuló hőmérséklet eloszlás az 2. peremfeltételek mellett [oC], 20µm-es légrés mellett Az 5.25. ábra a kialakult hőmérsékleteket mutatja a fűtőlap alatti légrés méretének függvényében, a két különböző peremfeltétel típus mellett.
75
5. Mikromechanikai szerkezetek kapcsolt modellezése
450
Hőmérséklet [Celsius]
400 350 300 250 200 150 100 50 0 20
40
60
80
100
Légrés [mikron]
PF1
PF2
5.25. ábra: Hőmérséklet a légrés méretének függvényében A hőterhelés miatt a szerkezetben hőfeszültségek keletkeznek. Ennek egyik oka a geometriai peremfeltételek, a másik pedig az, hogy a platina és a szilícium-nitrid hőtágulási együtthatója nagy mértékben különbözik. Mivel a platina teljesen be van ágyazva a szilícium-nitridbe, ezért a hőtágulása teljesen meg van gátolva. Az érzékelő középsíkjában kialakuló egyenértékű feszültség (HMH) az 5.26. ábrán látható. Az 5.27. ábra az elektromos hozzávezetésben kialakult egyenértékű feszültég eloszlását mutatja a fűtőlap alatti légrés méretének függvényében, a két különböző peremfeltétel típus mellett.
5.26. ábra: A Pt szálban kialakuló feszültség
76
5. Mikromechanikai szerkezetek kapcsolt modellezése
900
Egyenértékű feszültség [MPa]
800 700 600 500 400 300 200 100 0 20
40
60
80
100
Légrés [mikron] PF1
PF2
5.27. ábra: Egyenértékű feszültség (HMH) a légrés méretének függvényében Néhány eszköz el is készült, amelyek a felfűtés során valóban eltörtek a végeselem modell által is kiszámított hely környékén (5.28. ábra).
5.28. ábra: Repedések a megvalósított mintán A repedések elkerülése elsősorban a platina helyettesítésével érhető el, ezt a következő alfejezetben tárgyalom.
5. Mikromechanikai szerkezetek kapcsolt modellezése
77
5.2.2. Kéthidas gázérzékelő szerkezet A kéthidas gázérzékelőn az elektromos hozzávezetés anyagának megválasztását helyezem középpontba. Az előző pontban a peremfeltételek hatásának vizsgálatakor egy részletes modellt építettem fel, amellyel a levegőben történő hőterjedést is modellezni lehetett. Most azonban a szerkezetben kialakuló feszültségek alakulása az érdekes, ezért itt most csak az érzékelőn belüli hővezetést veszem figyelembe. A terhelés is más lesz, mint az előző esetben. A fűtőlap hőmérsékletét – a korábbi számítások alapján – 300 oC-nak veszem fel, a felfüggesztő hidak külső végei pedig szobahőmérsékletűek (20 oC). Azért, hogy a szilícium-nitridben kialakuló feszültségek valamelyest csökkenjenek, a felfüggesztő hidak egy átmeneti résszel kapcsolódnak a fűtőlaphoz. A modellezésnél szintén kihasználhatjuk a szerkezet szimmetriáját, azaz elegendő csupán a negyedét modellezni.
5.29. ábra: A gázérzékelő modell szerkezete Az 5.29. ábrán a gázérzékelő modellének szerkezetét láthatjuk, amely kismértékben eltér a valóságtól, hiszen a platina vagy a polikristályos szilícium elektromosan csatlakozik a szilíciumhoz. Azonban ez fölöslegesen bonyolította volna a modell elkészítését.
5.30. ábra: A gázérzékelő végeselem modellje
5. Mikromechanikai szerkezetek kapcsolt modellezése
78
Két anyagot vizsgálunk: a platinát és a polikristályos szilíciumot. Az 5.31. ábra a kialakult hőmérséklet eloszlást mutatja. A kialakult feszültségeket mutatja platina hozzávezetés esetén az 5.32. ábra, illetve polikristályos szilícium esetén az 5.33. ábra. Az ábrák alapján jól látható, hogy a polikristályos szilícium esetén kialakuló feszültségek sokkal kedvezőbbek, mint a platina esetében. Ennek alapvető oka az, hogy a polikristályos szilícium hőtágulási együtthatója (2,6⋅10-6 1/K) közel esik a szilícium-nitridéhez (2,11⋅10-6 1/K). Mivel így a 2. és 3. főfeszültségek is jóval kisebbek, ezért az egyenértékű feszültség eloszlása a szerkezetben jellegre megváltozik.
5.31. ábra: Hőmérséklet eloszlás a szerkezetben [oC]
5.32. ábra: Feszültség eloszlás a szerkezetben platina kontaktus esetén, a: egyenértékű feszültség, b: x irányú feszültség, c: y irányú feszültség, d: z irányú feszültség [kPa]
5. Mikromechanikai szerkezetek kapcsolt modellezése
79
5.33. ábra: Feszültség eloszlás a szerkezetben polikristályos szilícium kontaktus esetén, a: egyenértékű feszültség, b: x irányú feszültség, c: y irányú feszültség, d: z irányú feszültség [kPa] Ennél a kéthidas érzékelőnél a platina kontaktus esetén is kisebb feszültségek alakultak ki, mint a négyhidas érzékelő esetében, azonban ez is túl nagynak bizonyult a valóságos szerkezetnél (5.34. ábra). A polikristályos szilícium ezzel szemben megoldást jelenthet a nagy hőfeszültségek elkerülésére.
5.34. ábra: Repedések a kéthidas mintán
5. Mikromechanikai szerkezetek kapcsolt modellezése
80
A végleges szerkezeten mérést is végeztünk. Az 5.35. ábra a teljesítmény függvényében mutatja a felületen mért, illetve a végeselem modellel számított hőmérsékletet. A végeselem modell esetében csak a hővezetést vettük figyelembe, ez lehet a magyarázata a két görbe jelleg beli különbségének.
5.35. ábra: A megvalósított szerkezetn mért és számított hőmérséklet a teljesítmény függvényében
Új tudományos eredmény: Megmutattam, hogy mikromechanikai eszközök modellezése esetén a termikus hatások nem hanyagolhatóak el. Elkészítettem ezen eszközök kapcsolt termomechanikai modellét, melyekkel a valóságos viselkedés jobban leírható. A termikus mező hatását két különböző típusú kapacitív gyorsulásérzékelő szerkezeten a geometriai paramétereik függvényében mutattam meg. Gázérzékelő szerkezeteken mutattam be a peremfeltételek és az anyagjellemzők változtatásának hatását. A fejezet témájában megjelent publikációim: [P1], [P7], [P8], [P13], [P20], [P22] [P2], [P3], [P4], [P9], [P10], [P12], [P14], [P19], [P21]
6. Szilícium-nitrid rugalmassági moduluszának meghatározása
81
6. SZILÍCIUM-NITRID RUGALMASSÁGI MODULUSZÁNAK MEGHATÁROZÁSA Mikromechanikai rendszerek mechanikai modellezéséhez szükséges a felhasznált anyagok anyagjellemzőinek ismerete. Ezek között számos új anyag található. A mikromechanikai érzékelők legfontosabb alapanyagára, a szilíciumra, már elég megbízható anyagjellemzők találhatóak, sőt az esetleges adalékolások hatásáról is jelentek már meg közlemények [33], [35]. Más anyagokra azonban az irodalomban vagy nem, vagy egymástól jelentősen eltérő anyagjellemző értékek találhatóak az adott anyagra. A mikromechanikában széles körben alkalmazott nem sztöchiometrikus szilícium-nitrid tipikus példa erre. A szilícium-nitrid rendkívül vonzó tulajdonsága, hogy a szilícium és a nitrogén sztöchimetrikus arányának változtatásával a szilícium-nitrid elektromos, termikus és mechanikai tulajdonságai változtathatóak. Ez lehet az irodalomban található eltérések oka is. A legegyszerűbb mechanikai modellezéshez is szükséges anyagjellemző a rugalmassági modulusz. Ez az érték az irodalmi adatok alapján 100 és 400 GPa között található [77]. Azért, hogy az MTA Műszaki Fizikai és Anyagtudományi Intézetben készített érzékelők megfelelő pontossággal modellezhetőek lehessenek, a használt szilícium-nitrid rugalmassági moduluszát kellene meghatározni. A hagyományos mérési módszerek nem használhatók, hiszen szilícium-nitrid vékonyrétegről van szó. Egy elterjedt mikroméretekben használható megoldás a rezonancia módszer [33], [78]. A rezonancia módszer lényege, hogy a szerkezet sajátfrekvenciát mérjük meg, majd ebből számítjuk ki a rugalmassági moduluszt.
6.1.Mérés
6.1.1.Mérési elrendezés Kis méretű befogott tartóként modellezhető különböző geometriájú rudakat mértünk. A rudak egy chipen vannak kialakítva, ezt a chipet egy hangszóróra rögzítjük. A vizsgált rúd szabad végére lézer fényt (628 nm) bocsátottunk, amely onnan egy ernyőre verődik vissza. A rudakat a hangszóróval, hanggernerátorral gerjesztjük, ismert frekvenciával. A rezonancia jelensége a megfelelő gerjesztési frekvenciánál az ernyőn észlelhető, amikor a lézerfény nagy amplitúdójú rezgéseket végez. A frekvencia leolvasási pontossága ±25 Hz. A mérési elrendezés a 6.1.-6.2. ábrákon látható.
6. Szilícium-nitrid rugalmassági moduluszának meghatározása
6.1. ábra: A mérési elrendezés vázlata
6.2. ábra: A mérési elrendezés képe
82
6. Szilícium-nitrid rugalmassági moduluszának meghatározása
83
6.1.2. A minták A mérendő minták hagyományos vékonyréteg megmunkálási eljárással készültek. Szilíciumban gazdag szilícium-nitrid réteget (SiN1.05) hoztunk létre alacsony nyomású kémiai gőzöltető (LPCVD) eljárással, majd alkalikus anizotróp nedves marással alakítottuk ki a rudak alatti légrést. A 6.3. ábrán szilíciumból készült minta látható, mivel a szilíciumnitrid mintáról még nem készült kép.
6.3. ábra: Szilícium minta A minták különböző geometriával készültek. A rudak végén egy nagyobb, nyolcszögletű (a képen látható szilícium esetében négyzet alakú) platnit alakítottunk ki, aminek több feladata is van: • A számításokban az analitikus modell egy egyik végén befogott tartó, szabad végén koncentrált tömeggel. Ehhez a rúd szabad végén levő tömeg kialakítás előnyös • A lézer fény rögzített a mérés során. Ahhoz, hogy a fény végig a rúd végére essen nagyobb felület szükséges. • Az így kialakított nagyobb felületre további adalékok vihetők fel. Arany réteg felvitele tovább növeli a rúd végén levő tömeget. • Az arany a jó fényvisszaverő képessége miatt egyébként is szükséges a méréshez. Négy különböző geometriájú rudat vizsgáltunk, ezek adatait a 6.1. táblázat tartalmazza. Hossz Szélesség Vastagság Platni fő mérete Aranyréteg 100 30 0,7 200 0,1 300 30 0,7 200 0,1 100 60 0,7 200 0,1 300 60 0,7 200 0,1
6.1. táblázat: A rudak geometriája (µm)
6. Szilícium-nitrid rugalmassági moduluszának meghatározása
84
6.2. Számítások A mért eredmények alapján a sajátfrekvenciából számítással lehet a rugalmassági moduluszt meghatározni. Ehhez két különböző modellt használtam. Az első egy analitikus rud modell, míg a másik egy végeselem modell.
6.2.1. Analitikus modell Az analitikus rúdmodell esetén a rúd maga csak a rugó szerepét látja el, tömegét elhanyagoltam, a rúd végén levő koncentrált tömeg nagysága a sűrűség adatokból és a geometriából számolható ki.
6.4. ábra: Analitikus rúdmodell A rúd, mint rugó merevségét a rúdvég elmozdulása alapján lehet meghatározni: s=
3IE , l3 (6.1)
ahol s – a rugó merevsége, I – a rúd keresztmetszetének másodrendű nyomatéka, E – a rugalmassági modulusz, l – a rúd hossza. Egy szabadságfokú modellről van szó. A sajátkörfrekvencia és abból a sajátfrekvencia könnyen kiszámítható: α=
s , m
(6.2) f=
α , 2π (6.3)
6. Szilícium-nitrid rugalmassági moduluszának meghatározása
85
ahol α – a sajátkörfrekvencia, f – a sajátfrekvencia, m – a tömeg.
Ezek alapján a rugalmassági modulusz meghatározására szolgáló képlet: E=
4f 2 ⋅ m ⋅ l 3 ⋅ π 2 . 3I (6.4)
6.2.2. Végeselem modell A végeselem modell elkészítéséhez a COSMOS/M 2.0 programot használtam. A modell a 6.5. ábrán látható.
6.5. ábra: Végeselem modell
6.2.3. Eredmények A rugalmassági moduluszt alapvetően a végeselem modell segítségével számítottam ki, a rúdmodellt becslésre használtam a végeselem modell bemenő paraméteréhez. A 6.2. táblázatban szereplő eredmények az átlagos értékek az adott geometriára.
86
6. Szilícium-nitrid rugalmassági moduluszának meghatározása
Hossz
Szélesség Vastagság
Tömeg [kg]
Frekvencia [Hz]
Rugalmassági modulusz [GPa] Analitikus VEM 202 202
100
30
0,7
1,31 ⋅ 10 −10
3550
300
30
0,7
1,31 ⋅ 10 −10
1250
201
204
100
60
0,7
1,31 ⋅ 10 −10
4750
181
183
300
60
0,7
1,31 ⋅ 10 −10
1750
197
202
6.2. táblázat: A szilícium-nitrid rugalmassági modulusza A fenti eredményeket referencia eredményekkel is összehasonlítottam. A szilícium rugalmassági modulusza ismert, ezért szilícium rudakra is elvégeztük a fenti méréssorozatot. Ezáltal a fenti módszer megbízhatósága ellenőrízhető. Az N-típusú szilícium rudakat porózus szilícium felhasználásával készítettük. A 6.6. ábrán látható, hogy a mért és számított eredmények jó egyezést mutatnak.
4,00E+04 3,50E+04 Frekvencia [Hz]
3,00E+04 2,50E+04 2,00E+04 1,50E+04 1,00E+04 5,00E+03 0,00E+00 0,00E+00
2,00E-04
4,00E-04
6,00E-04
8,00E-04
1,00E-03
Rúd hossza [m] Analitikus számítás
Mérés
Végeselem számítás
6.6. ábra: A referencia mérés eredményei Az ábrán két sorozat látható, az alsó a 30 µm, a felső a 60 µm széles rúdhosszhoz tartozik. A fentiek alapján megállapítható, hogy a módszer használható a rugalmassági modulusz meghatározására. A szilícium-nitrid (SiN1.05) rugalmassági modulusza 180-210 GPa tartományba esik. Ez az intervallum elegendően szűk ahhoz, hogy szilícium-nitridet tartalmazó mikromechanikai eszközök mechanikai modellezéshez felhasználható legyen.
6. Szilícium-nitrid rugalmassági moduluszának meghatározása
87
Új tudományos eredmény: A szilicium-nitrid mechanikai modellezéséhez feltétlenül szükséges rugalmassági modulusz meghatározását frekvencia mérés után analitikus és végeselem módszer segítségével végeztem el. A mérés során kis méretű befogott rúdként modellezhető rudakat hangfrekvenciás tartományban gerjesztettünk. A rudak rezonancia frekvenciája alapján az analitikus és végeselem modellekkel kiszámítottam a szilícium-nitrid rugalmassági moduluszát. A fejezet témájában megjelent publikációim: [P17]
Összefoglalás
88
ÖSSZEFOGLALÁS Munkámban mikromechanikai eszközök modellezése során felmerülő kapcsolt feladatok megoldásával foglalkoztam. A mikromechanikai szerkezetekben leggyakrabban mechanikai, elektromos és termikus hatások fordulnak elő. Levezettem a kapcsolt piezo-termo-elasztikus anyag- és mezőegyenleteket egy módosított hiperbolikus hővezetési egyenlet felhasználásával. A Fourier-törvény alkalmazásával a kapcsolt piezo-termo-elasztikus egyenletek már ismertek. A hiperbolikus hővezetési egyenlet felhasználása azért fontos, mert a hőterjedéssel kapcsolatos jelenségek így nem végtelen nagy sebességgel terjednek, és ennek a mikronos méret tartományban nagy jelentősége van. Galjorkin módszerének felhasználásával előállítottam a kapcsolt piezo-termo-elasztikus egyenletek végeselemes diszkretizációját: • a klasszikus, Fourier-féle hővezetési egyenlet alkalmazásával. Az irodalomban már találhatóak ilyen végeselemes levezetések, azonban ezek egymástól eltérőek és egyiket sem a Galjorkin módszerrel vezeték le. • a módosított, hiperbolikus hővezetési egyenlet felhasználásával. A piezo-termo-elasztikus egyenletek végeselemes diszkretizációjával előállított végeselemes mátrix-differenciálegyenletre vonatkozóan az irodalomban már található egykét példa, azonban ezek különbözőek. Az általam Galjorkin módszerével előállított egyenlet az utóbbi két [20] és [21] cikk eredményeit támasztja alá. A fenti egyenletek alapján végeselem programot készítettem. A kapott eredményeket ismertettem, összevetettem őket más végeselem programok eredményeivel. A kereskedelmi végeselem programok nem tudnak három kapcsolt mezőt kezelni (legfeljebb páronként), így az összevetésből látható eltérés indokolja a piezo-termo-elasztikus egyenletek használatát. Mintafeladatot mutatok be, amelyen jól látható az eltérés a mechanikai, termomechanikai, piezoelektromos és piezo-termo-elasztikus analízisek eredményei között. A programban a mátrix egyenletet egy lépésben oldom meg (ún. direkt módon). Megmutattam, hogy mikromechanikai eszközök modellezése esetén a termikus hatások nem hanyagolhatóak el. Elkészítettem ezen eszközök kapcsolt termomechanikai modellét, melyekkel a valóságos viselkedés jobban leírható. Bemutattam a termikus mező hatását két különböző típusú kapacitív gyorsulásérzékelő szerkezeten a geometriai paramétereik függvényében. Gázérzékelő szerkezeteken mutattam be a peremfeltételek és az anyagjellemzők változtatásának hatását. A szilícium-nitrid a mikromechanikában gyakran alkalmazott anyag, rossz elektromos és hővezetése miatt. Mechanikai tulajdonságai nem ismertek kellő pontossággal. A mechanikai modellezéshez feltétlenül szükséges rugalmassági modulusz meghatározását analitikus és végeselem módszer segítségével végeztem el.
Summary
89
SUMMARY The topic of this work is modelling coupled mechanical problems, which occurs in MEMS (Micro-Electro-Mechanical Systems). In the most cases one has to deal with mechanical, thermal and electric effects. In the first parts in my work the derivation and finite element discretization of piezothermoelastic equations are presented. Numerical calculations were carried out on modelling examples with an own finite element code. The next part showed the thermomechanical analysis of real micromechanical structures. Finally, elasticity modulus measurement of silicon-nitride is presented. I derived the coupled piezothermoelastic equations using a modified hyperbolic equation of heat conduction (suggested by CATTANEO and VERNOTTE). The hyperbolic equation of heat conduction is more adequate in micro scale. The equations are valid for homogenous, linear elastic, anisotropic material. I discretized the piezothermoelastic equations using Galerkin’s method with • the classical Fourier’s law • the modified hyperbolic equation of heat conduction There are a few publications dealing with finite element discretization of piezothermoelastic equations, but the presented equations are different form each other. My results confirm the equations of GÖRNANDT et al. and BOURNATA et al. Both were derived using variation principles. The hyperbolic equation of heat conduction was not considered in any publication. I developed an own finite element code in Fortran to solve coupled piezothermoelastic problems. I used direct solution strategy to solve the equations. Sequential strategy can be applied too, but the direct strategy allows to take more difficult material laws into account in the future. The disadvantage of the direct method is that the stiffness matrix is unsymmetrical. I presented results of thermomechaical analysis for real micromechanical structures. The effects of geometrical parameters on the behaviour of structure are showed on a capacitive type accelerometer. The effect of boundary conditions and material types are presented on an isolated, suspended gas sensor. The elasticity modulus measurement of deposited silicon-nitride layer was described. Silicon-nitride is a widely used structural material in MEMS applications because of its electrical and thermal properties. But its mechanical properties are given in a wide range in the literature. The Young’s modulus can be determined from the first natural frequency of thin cantilever beams formed from silicon-nitride layers. The elasticity modulus, which is the most important mechanical parameter for structural modelling, was calculated both from analytical and finite element models.
90
Tézisek
TÉZISEK 1.
2.
3.
4.
5.
Levezettem a kapcsolt piezo-termo-elasztikus anyag- és mezőegyenleteket egy módosított (CATTANEO és VERNOTTE által javasolt), hiperbolikus hővezetési egyenlet felhasználásával. A Fourier-törvény alkalmazásával a kapcsolt piezotermo-elasztikus egyenletek már ismertek. A hiperbolikus hővezetési egyenlet felhasználása azért fontos, mert a hőterjedéssel kapcsolatos jelenségek így nem végtelen nagy sebességgel terjednek, és ennek a mikronos mérettartományban jelentősége van. Az egyenletek lineárisan rugalmas, homogén, anizotróp esetben érvényesek, kis alakváltozás mellett. A mezőegyenletek változói az elmozdulás, a hőmérséklet-különbség és az elektromos potenciál. Ezek lehetnek a végeselemes egyenletek változói is. Galjorkin módszerének felhasználásával előállítottam a kapcsolt piezo-termoelasztikus egyenletek végeselemes diszkretizációját: • a klasszikus, Fourier-féle hővezetési egyenlet alkalmazásával. • a módosított, hiperbolikus hővezetési egyenlet felhasználásával. Mindkét esetre levezettem a diszkretizált mátrix differenciál egyenletet, amely alapján végeselemes algoritmus készíthető. A 2. tézisben levezetett egyenletek alapján végeselem algoritmust és programot készítettem. A programban a mátrix egyenletet - a korábbi publikációkban közöltekkel ellentétben, ahol a termikus részt leválasztva két lépésben (ún. szekvenciálisan) oldják meg az egyenleteket - egy lépésben (ún. direkt módon) oldom meg. Megmutattam, hogy mikromechanikai eszközök modellezése esetén a termikus hatások nem hanyagolhatóak el. Elkészítettem ezen eszközök kapcsolt termomechanikai modelljét, melyekkel a valóságos viselkedés jobban leírható. A termikus mező hatását két különböző típusú kapacitív gyorsulásérzékelő szerkezeten a geometriai paramétereik függvényében mutattam meg. Gázérzékelő szerkezeteken mutattam be a peremfeltételek és az anyagjellemzők változtatásának hatását. A szilicium-nitrid mechanikai modellezéséhez feltétlenül szükséges rugalmassági modulusz meghatározását frekvencia mérés után analitikus és végeselem módszer segítségével végeztem el. A mérés során kis méretű befogott rúdként modellezhető rudakat hangfrekvenciás tartományban gerjesztettünk. A rudak rezonancia frekvenciája alapján az analitikus és végeselem modellekkel kiszámítottam a szilícium-nitrid rugalmassági moduluszát.
Irodalom
91
IRODALOM [1] G. L. Pearson, W. T. Read Jr., and W. L. Feldmann: „Deformation and Fracture of Small Silicon Crystals” Acta Metallurgica, Vol. 5, 181-191. (1957) [2] R. D. Mindlin: „On the equations of motion of piezoelectric crystals”, Problems of Continuum Mechanics, Muskhelishvili, N. I., 70th birthday volume, pp. 282-290. (1961) [3] J. L. Nowinski: „Theory of thermoelasticity with applications”, Sijthoff & Noordhoff International Publishers B.V., Alphen aan den Rijn, The Netherlands (1978) [4] G. A. Altay, M. C. Dökmeci: „Some comments on the higher order theories of piezoelectric, piezothermoelastic and thermopiezoelectric rods and shells”, International Journal of Solids and Structures 40, 4699-4706. (2003) [5] R. D. Mindlin: „Equations of High Frequency Vibrations of Thermopiezoelectric Crystal Plates”, International Journal of Solids and Structures, Vol. 10, pp. 625-637. (1974) [6] G. P. Dube, S. Kapuria, P. C. Dumir “Exact Piezothermoelastic Solution of SimplySupported Orthotropic Flat Panel In Cylindrical Bending”, Int. J. Mech. Sci., Vol. 38, No. 11, 1161-1177. (1996) [7] H.-J. Lee, D. A. Saravanos: „A Mixed Multi-Field Finite Element Formulation for Thermopiezoelectric Composite Shells”, SPIE Conference on Mathematics and Control in Smart Structures, Newport Beach, California, SPIE Vol. 3667, pp. 449-460. (1999) [8] F. Ashida, T. R. Tauchert, „A general plane-stress solution in cylindrical coordinates for a piezothermoelastic plate”, International Journal of Solids and Structures 38, 4969-4985. (2001) [9] M. Stam, G. Carman: „Electrothermoelastic Behavior of Piezoelectric Coupled Cylinders”, American Institute of Aeronautics and Astronautics Journal, Vol. 34, No. 8, 1612-1618. (1996) [10] C.-Q. Chen, Y.-P. Shen: „Piezothermoelasticity Analysis for a Circular Cylindrical Shell Under the State of Axisymmetric Deformation”, International Journal of Engineering Science, Vol. 34, No. 14, pp. 1585-1600. (1996)
Irodalom
92
[11] Q. Qin: „Thermoelectroelastic analysis of cracks in piezoelectric half-plane by BEM”, Computational Mechanics 23, 353-360. (1999) [12] Q.-H. Qin, Y.-W. Mai: „Thermoelectroelastic Green’s function and its application for bimaterial of piezoelectric materials”, Archive of Applied Mechanics 68, 433-444. (1998) [13] Q.-H. Qin: „Green function and its application for a piezoelectric plate with various openings”, Archive of Applied Mechanics 69, 133-144. (1999) [14] Q.-H. Qin: „Green’s function for thermopiezoelectric plates with holes of various shapes”, Archive of Applied Mechanics 69, 406-418. (1999) [15] A. K. Belyaev: „Thermodynamic rationale for heat conduction equation and dynamic boundary value problem for piezothermoelastic materials”, SPIE Conference on Mathematics and Control in Smart Structures, Newport Beach, California, SPIE Vol. 3667, pp. 702-710. (1999) [16] J.-Q. Tarn: „A state space formalism for piezothermoelasticity”, International Journal of Solids and Structures 39, 5173-5184. (2002) [17] R. V. N. Melnik: „Computational analysis of coupled physical fields in piezothermoelastic media”, Computer Physics Communications 142, 231-237. (2001) [18] R.-F. Fung, S.-C. Chao, Y.-S. Kung: „Piezothermoelastic analysis of an optical beam deflector”, Sensors and Actuators A, 87, 179-187. (2001) [19] T. S. Koko, M. J. Smith, I. R. Orisamolu: „Modelling Stresses in Piezoelectric Smart Structures under Combined Thermal and Mechanical Excitations”, SPIE Conference on Mathematics and Control in Smart Structures, Newport Beach, California, SPIE Vol. 3667, pp. 474-485. (1999) [20] K. Bouranta, G. A. Malegiannakis, B. Kröplin: „Thermo-electro.mechanical coupling problems solved by FE-formulation”, Engineering Computations, Vol. 15 No. 6, pp. 804-828. (1998) [21] A. Görnandt, U. Gabbert: „Finite element analysis of thermopiezoelectric smart structures”, Acta Mechanica 154, 129-140 (2002) [22] Balla M., „A hőrugalmasságtan kapcsolt feladatainak analitikus és végeselemes megoldásai”, kandidátusi értekezés, Budapest (1994)
Irodalom
93
[23] ANSYS, ANSYS, Inc., Southpointe, 275 Technology Drive, Canonsburg, PA 15317, USA [24] J. Funk, J. G. Korvink, M. Bächtold, J. Bühler, H. Baltes: „Coupled 3D Thermoelectro-mechanical Simulations of Microactuators”, Proceedings of the IEEE Micro Electro Mechanical Systems (MEMS), p. 133-138. (1996) [25] G.-S. Chen, M.-S. Ju, Y.K. Fang: „Effects of monolithic silicon postulated as an isotropic material on design of microstructures”, Sensors and Actuators A, 86, 108-114. (2000) [26] E. Mazza, J. Dual: „Mechanical behavior of a µm-sized single crystal silicon structure with sharp notches”, Journal of the Mechanics and Physics of Solids 47, 17951821. (1999) [27] L. Wei, M. Vaudin, C. S. Hwang, G. White: „Heat conduction in silicon thin films: Effect of microstructure”, J. Mater. Res., Vol. 10, No. 8, 1889-1896. (1995) [28] W. Fang, H.-C. Tsai, C.-Y. Lo: „Determining thermal expansion coefficients of thin films using micromachined cantilevers”, Sensor and Actuators A, 77, 21-27. (1999) [29] W. Fang, C.-Y. Lo: „On the thermal expansion coefficients on thin films”, Sensor and Actuators A, 84, 310-314. (2000) [30] V. Ziebart, O. Paul, H. Baltes: „Extraction of the coefficient of thermal expansion of thin films from buckled membranes”, Mat. Res. Soc. Symp. Proc. Vol. 546, 103-109. (1999) [31] K. E. Petersen: „Silicon as a Mechanical Material”, Proceedings of the IEEE, vol. 70, No. 5., 420-457. (1982) [32] C. S. Pan, W. Hsu: „A Microstructure for in situ Determination of Residual Strain”, Journal of Microelectromechanical Systems, vol. 8, No. 2 (1999) [33] S. Lee, C. Cho, J. Kim, S. Park, S. Yi, J. Kim, D. D. Cho: „The effects of postdeposition processes on polysilicon Young’s modulus”, Journal of Micromechanics and Microengineering 8 330-337. (1998) [34] W-H. Chu, M. Mehregany: „A Study of Residual Stress Distribution Through the Thickness of p+ Silicon Films”, IEEE Transactions on Electron Devices, vol. 40. NO. 7, 1245-1250. (1993)
Irodalom
94
[35] X. Ding, W. H. Ko, J. M. Mansour: „Residual Stress and Mechanical Properties of Boron-doped p+-Silicon Films” , Sensors and Actuators A, 21-23, 866-871. (1990) [36] F. Rudolf: „A Micromechanical Capacitive Accelerometer with a Two-point Internal-mass Suspension” Sensors and Actuators A, 4, 191-198. (1983) [37] H. V. Allen, S. C. Terry, and D. W. de Bruin: „Accelerometer Systems with Selftestable Features” Sensors and Actuators A, 20, 153-161. (1989) [38] R. Puers, A. Montano, S. Reyntjens: „Design of a new miniaturized capacitive triaxial accelerometer” Proceedings of the EUROSENSORS XI, Warsaw, Poland, September 21-24, 755-758. (1997) [39] R. Puers, S. Reyntjens: „The characterization of a miniature silicon micromachined capacitive accelerometer”, Journal of Micromechanics and Microengineering 8, 127-133. (1998) [40] V. Josselin, P. Touboul, R. Kielbasa: „Capacitive detection scheme for space accelerometers applications”, Sensors and Actuators A, 78, 92-98. (1999) [41] W. Qu, C. Wenzel, G. Gerlach: „Fabrication of a 3D differential-capacitive acceleration sensor by UV-LIGA”, Sensors and Actuators A, 77, 14-20. (1999) [42] D. Haronian: „A low-cost micromechanical accelerometer with integrated solidstate sensor”, Sensors and Actuators A, 84, 149-155. (2000) [43] O. Bochobza-Degani, D. J. Seter, E. Socher, Y. Nemirovsky: „A novel micromachined vibrating rate-gyroscope with optical sensing and electrostatic actuation”, Sensors and Actuators A, 83, 54-60. (2000) [44] D. L. Howard, S. D. Collins, R. L. Smith: „A single-fringe etalon silicon pressure transducer”, Sensors and Actuators A, 86, 21-25. (2000) [45] D. De Bruyker, R. Puers: „Thermostatic control for temperature compensation of a silicon pressure sensor”, Sensors and Actuators A, 82, 120-127. (2000) [46] L. M. Roylance, J. B. Angell: „A Batch-Fabricated Silicon Accelerometer” IEEE Trans. on El. Dev. Vol. Ed-26. No. 12. 1911-1917. (1979) [47] H. Seidel, and L. Csepregi: „Design Optimization for Cantilever-type Accelerometers” Sensors and Actuators A, 6 81-92. (1984)
Irodalom
95
[48] S. Shen, J. Chen, and M. Bao: „Analysis on twin-mass structure for a piezoresistive accelerometer” Sensors and Actuators A, 34 101-107. (1992) [49] H. E. A. Elgamel: „A simple and efficient technique for the simulation of capacitive pressure transducers”, Sensors and Actuators A, 77, 183-186. (1999) [50] L. Cao, T. S. Kim, S. C. Mantell, D. L. Polla: „Simulation and fabrication of piezoresistive membrane type MEMS strain sensors”, Sensors and Actuators A, 80, 273279. (2000) [51] H. Crazzolara, G. Flach, W. von Münch: „Piezoresistive accelerometer with overload protection and low cross-sensitivity”, Sensors and Actuators A, 39, 201-207. (1993) [52] K. Yoshida, Y. Matsumoto, M. Ishida, K. Okada: „High-sensitive three axis SOI capacitive accelerometer using dicing method”, Tech. Dig. of the 16th Sensor Symposium, pp.25-28. (1998) [53] O. Lüdtke, V. Biefeld, A. Buhrdorf, J. Binder: „Laterally driven accelerometer fabricated in single crystalline silicon”, Sensors and Actuators A, 82, 149-154. (2002) [54] Z. Li, Z. Yang, Z Xiao, Y. Hao, T. Li, G. Wu, Y. Wang: „A bulk micromachined vibratory lateral gyroscope fabricated with wafer bonding and deep trench etching”, Sensors and Actuators A, 83, 24-29. (2000) [55] A. Chouaf, Ch. Malhaire, M. Le Berre, M. Dupeux, F. Pourroy, D. Barbier: „Stress analysis at singular points of micromachined silicon membranes”, Sensors and Actuators A, 84, 109-115. (2000) [56] Z. Xiao, M. Chen, G. Wu, C. Zhao, D. Zhang, Y. Hao, G. Zhang, Z. Li: „Silicon micro-accelerometer with mg resolution, high linearity and large frequency bandwith fabricated with two mask bulk process”, Sensors and Actuators A, 77, 113-119. (1999) [57] J. A. Plaza, J. Esteve, C. Cané: „Twin-mass accelerometer optimization to reduce the package stresses”, Sensors and Actuators A, 80, 199-207. (2000) [58] G. Blasquez, X. Chauffleur, P. Pons, C. Douziech, P. Favaro, Ph. Menini: „Intrinsic thermal behaviour of capacitive pressure sensors: mechanisms and minimisation”, Sensors and Actuators A, 85, 65-69. (2000) [59] R. Puers, and D. Lapadatu: „Electrostatic forces and their effects on capacitive mechanical sensors” Sensors and Actuators A, 56 203-210. (1996)
Irodalom
96
[60] R. Luck, E. I. Agba: „On the design of piezoelectric sensors and actuators”, ISA Transactions 37, 65-72. (1998) [61] J. A. Hauch, D. Holland, M. P. Marder, H. L. Swinney: „Dynamic Fracture in Single Crystal Silicon”, Physical Review Letters Vol. 82, No. 19, 3823-3826. (1999) [62] W. N. Sharpe Jr., K. M. Jackson, K. J. Hemker, Z. Xie: „Effect of Specimen Size on Young’s Modulus and Fracture Strength of Polysilicon”, Journal of Microelectromechanical Systems, Vol. 10, No. 3, 317-326. (2001) [63] N. Tayebi, A. K. Tayebi: „Numerical Fracture Analysis of MEMS devices”, International Conference on Modeling, Simulation of Microsystems MSM99. Puerto Rico, April 99. (1999) [64] C. Muhlstein, S. Brown: „Reliability and Fatigue Testing of MEMS”, NDF/AFOSR/ASME Workshop, Tribology Issues and Opportunities in MEMS. Nov. 9-11 (1997) [65] U. M. Mescheder, W. Kronast, N. Naychuk: „Reliability Investigations in micromechanical Devices”, The 16th European Conference on Solid-State Transducers, September 15-18, Prague, Czech Republic (2002) [66] I. M. Smith, D. V. Griffiths „Programming the Finite Element Method”, John Wiley & Sons Ltd., Chichester, England (1998). [67] COSMOS, Structural Research & Analysis Corporation (SRAC), 12121 Wilshire Blvd., Suite 700, Los Angeles, CA 90025 [68] C. Cattaneo: „A Form of Heat Equation which Eliminates the Parados of Instantaneous Propagation”, CR Acad. Sci. 247, 431-433 (1958). [69] P. Vernotte: „Paradoxes in the Continous Theory of the Heat Equation”, CR Acad. Sci. 246, 3154-3155 (1958). [70] M. Chester: „Second Sound in Solids”, Physical Review 131, 2013-2015 (1963). [71] P. H. Francis: „Thermomechanical Effects in Elastic Wave Propagation”, A Survey. Journal Sound Vibration 21, 181-192 (1972).
Irodalom
97
[72] D. Cahill, K. Goodson, A. Majumdar: „Thermometry and Thermal Transport in micro/nano Scale Solid-state Devices and Structures”, Journal Heat Transfer, ASME, Vol. 124, 223-241 (2002) [73] Gy. Gróf, L.I. Kiss: „Heat Conduction Models”, Proceedings of Fourth Conference on Mechanical Engineering, Gépészet 2004, Budapest, pp.353-356 (2004) [74] E. Zanchini: „Hyperbolic-heat-conduction Theories and Nondecreasing Entropy”, Physical Review B, Volume 60, Number 2, pp. 991-997, (1999) [75] A. Szekeres: „The Second Sound Phenomenon: Pro and Contra”, Periodica Polytechnica 48, No. 1, pp.83-87, (2004) [76] „Sensors” Edited by W. Göpel, J. Hesse, J.N. Zemel Vol. 7 „Mechanical Sensors” Edited by H. H. Bau, N. F. de Rooij, and B. Kloeck, VHC Weinheim, (1994) [77] http://www.memsnet.org/material/ [78] L. Kiesewetter and J.-M. Zhang D. Houdeau and A. Steckenborn: „Determination of Young's moduli of micromechanical thin films using the resonance method”, Sensors and Actuators A, 35, pp.153-159. (1992) [79] „Semiconductor sensors” Edited by S. M. Sze, Wiley, New York, 1994
Publikációk
98
PUBLIKÁCIÓK Lektorált folyóirat cikk [P1] Á. Kovács, Zs. Vízváry: "Structural parameter sensitivity analysis of cantileverand bridge- type accelerometers" Sensors and Actuators A, Vol. 89/3, 2001, pp.197-205 [P2] Zs. Vízváry, P. Fürjes, I. Bársony: "Thermomechanical Analysis of Hotplates by FEM" Microelectronics Journal 32, 2001, pp. 833-837 [P3] P. Fürjes, Zs. Vízváry, M. Ádám, A. Morrissey, Cs. Dücsõ, I. Bársony: "Thermal investigation of micro-filament heaters", Sensors and Actuators A, Vol. 99/1-2, 2002, pp.98-103 [P4] P. Fürjes, Zs. Vízváry, M. Ádám, I. Bársony, A. Morrissey, Cs. Dücsõ : "Materials and processing for realization of micro-hotplates operated at elevated temperature", Journal of Micromechanics and Microengineering, Vol. 12, No. 4 pp.425-429 (2002) [P5] I. E. Lukács, Zs. Vízváry, P. Fürjes, F. Riesz, Cs. Dücsõ and I. Bársony: "Determination of Deformation Induced by Thin Film Residual Stress in Structures of Millimeter Size", Advanced Engineering Materials, Vol. 4 – No. 8 pp.625-627 (2002) [P6] Vízváry Zs.: “Piezo-termo-elasztikus feladatok végeselemes modellezése”, Elektronika Technológia Mikrotechnika, 2004/2-3 (elfogadva)
Elektronikus publikáció [P7] Á. Kovács, Zs. Vízváry: "Analytical calculation of sensitivity for capacitive pressure transducers" PAMM, Vol. 1, Issue 1, pp. 135-136, (2002). Lektorált konferencia cikk [P8] Zs. Vízváry: "Modelling of silicon based accelerometer structures" Proceedings of Second Conference on Mechanical Engineering, Gépészet 2000, Budapest, 2000, pp.789793 [P9] Zs. Vízváry, P. Fürjes, I. Bársony: "Three-Dimensional Finite Element Model for Thermomechanical Analysis of Hotplates" Proceedings of the 6th International Workshop on THERMal INvestigations of ICs and Systems (Therminic), Budapest, 2000, pp. 253257 [P10] P. Fürjes, Zs. Vízváry, M. Rácz, I. Bársony: "Temperature measurement in micro-filament heater" Proceedings of the 6th International Workshop on THERMal INvestigations of ICs and Systems (Therminic), Budapest, 2000, pp. 262-265
Publikációk
99
[P11] I.E. Lukács, Zs. Vízváry, P. Fürjes, F. Riesz, Cs. Dücsõ, I. Bársony: "Determination of Deformation Induced by Thin Film Residual Stress in Structures of Millimetre Size" Proceedings of E-MRS 2001, Strasbourg, 2001, pp.O-8 (O/P03) [P12] P. Fürjes, Zs. Vízváry, M. Ádám, M. Rácz, I. Bársony: "Thermal investigation of Micro-filament Heaters" Proceedings of E-MRS 2001, Strasbourg, 2001, pp.J-14 (J-IV.10) [P13] Zs. Vízváry, Á. Kovács: "Design optimization for cantilever- and bridge-type acceleration microsensors using the Finite Element Method" European Conference on Computational Mechanics, Krakkó, 2001 (CD-ROM) [P14] Zs. Vízváry, P. Fürjes: "Thermomechanical Investigation of a Suspended Microhotplate" Proceedings of Third Conference on Mechanical Engineering, Gépészet 2002, Budapest, 2002, pp.302-306 [P15] Zs. Vízváry: "Finite Element Formulation of Piezothermoelastic Equations" Proceedings of the Fifth World Congress on Computational Mechanics (WCCM V), July 7-12, 2002, Vienna, Austria, Editors: Mang, H.A.; Rammerstorfer, F.G.; Eberhardsteiner, J., Publisher: Vienna University of Technology, Austria, ISBN 39501554-0-6, http://wccm.tuwien.ac.at [P16] Zs. Vízváry, P. Fürjes, M. Ádám, Cs. Dücsõ and I. Bársony:"Mechanical Modelling of an Integrable 3D Force Sensor by Silicon Micromachining", Proceedings of MME'02, Micromechanics Europe, Sinaia, 2002, pp.165-168 [P17] Zs. Vízváry, J. Volk, Cs. Dücsõ: "Elasticity Modulus Measurement Of Silicon-Nitride" Proceedings of Fourth Conference on Mechanical Engineering, Gépészet 2004, Budapest, 2004, pp.215-219
Egyéb konferencia cikk [P18] Zs. Vízváry, M. Ádám, I. Bársony: “Geometric Considerations In Micromachined Capacitive Pressure Sensors” International Workshop on MicroDevices, Budapest, 1999 [P19] P. Fürjes, Zs. Vízváry, M. Ádám, Cs. Dücsõ, A. Tóth, I. Bársony: "Processing and characterisation of integrable microhotplates for gas sensing applications" Book of Abstract, First Conference on Microelectronics, Microsystems, Nanotechnology, MMN 2000, Athén, 2000
Publikációk
100
Előadás [P20] Á. Kovács, Zs. Vízváry: "Analytical calculation of sensitivity for capacitive pressure transducers" GAMM Conference'01, Zürich, Switzerland, February, 12-15, 2001. [P21] Zs. Vízváry: "Mechanical Analysis and Optimisation of a Microhotplate" FinnoUgric International Conference of Mechanics, 27 May - 2 June, 2001. [P22] Zs. Vízváry, Á. Kovács: "Thermomechanical modelling of silicon based accelerometers" 3rd Korean-Hungarian Symposium on "New methods in structural engineering", Budapest, 2001 [P23] Zs. Vízváry: "Piezo-termo-elasztikus egyenletek végeslemes modellezése" IX. Magyar Mechanikai Konferencia, Miskolc 2003. augusztus 27-29.
Köszönetnyilvánítás
101
KÖSZÖNETNYILVÁNÍTÁS Munkám során sokan segítettek. Ezúton szeretném köszönetemet kifejezni témavezetőmnek, Kovács Ádámnak, az MTA Műszaki Fizikai és Anyagtudományi Kutató Intézetének, elsősorban Ádám Antalnének, Bársony Istvánnak és Dücső Csabának. Fürjes Péternek és Szabó Zsoltnak. A Tateyama Hungary Kft.-nek, Morita Tsuneo-nak. Az FHFurtwangennek, Kovács Andrásnak és Ulrich Meschedernek.