PORAZDELITEV TRDNOSTI KERAMIČNIH MATERIALOV Milan Ambrožič Institut "Jožef Stefan", Jamova 39, 1000 Ljubljana POVZETEK Meritve nekaterih mehanskih lastnosti keramike dajejo relativno {iroko porazdelitev vrednosti teh veli~in za serijo enakih vzorcev. Zgled je meritev upogibne trdnosti z enoosnim 3- ali 4-to~kovnim upogibnim preskusom. Porazdelitev trdnosti se navadno dobro opi{e z Weibullovo porazdelitveno funkcijo, in ni rezultat nenatan~nosti merjenja, temve~ naklju~ne porazdelitve napak v snovi. V prispevku je opisana simulacija Monte Carlo Weibullove porazdelitve. Strength distribution of ceramic materials ABSTRACT Measurements of some mechanical properties of ceramics result in relatively broad distribution for the series of equal samples. An example is the measurement of the bending strength with a 3- or 4-point bending test. The strength distribution can usually be well described by a Weibull distribution function. It is not the conseqence of the measurement uncertainty, but results from the random distribution of flaws in the material. The Monte Carlo simulation of the Weibull dristribution is presented in this paper. 1 UVOD Weibullova porazdelitev je znana že od sredine prejšnjega stoletja in je osnovana na principu "najšibkejšega člena", to je, da se material oz. izdelek zlomi, ko popusti njegov najšibkejši del (1'2). Weibullova porazdelitvena funkcija, ki v najosnovnejši obliki vsebuje dva prosta parametra, Weibullov modul in karakteristični parameter, je bila neštetokrat eksperimentalno potrjena na različnih področjih, od mikroelektronskih komponent do gradbenih materialov (2). Najuporabnejša je za opis porazdelitve dveh veličin: 1) trajnostnega časa enakih izdelkov iz serijske proizvodnje (časa od izdelave do odpovedi ali okvare izdelka) pri vsaj približno podobnih razmerah delovanja, 2) sile, pri kateri se material/izdelek zlomi (oz. iz sile in geometrije izračunane mehanske napetosti v snovi). Zaradi velikega pomena statistične zanesljivosti ocenjevanja vedenja materialov/izdelkov na osnovi dane porazdelitvene funkcije, ki jo praviloma izračunamo na osnovi relativno majhnega števila preskusnih vzorcev (vsak pri preskusu polomljeni vzorec pomeni strošek), je bilo narejenih veliko teoretičnih raziskav Weibullove porazdelitve (3-9). Priljubljena metoda je simulacija Monte Carlo, pri kateri si pomagamo z računalniškim generatorjem naključnih števil med 0 in 1. Vendar pa je bila doslej velika večina tovrstnih raziskav usmerjena k simulaciji same Weibullove porazdelitve, to je k njeni pravilni matematični interpretaciji, ne glede na njeno fizikalno ozadje. Cilj takšnih raziskav je doseči čim bolj zanesljivo napoved Weibullovih parametrov na osnovi majhnega števila meritev na vzorcih. Le majhno število raziskav je bilo namenjenih iskanju zveze med porazdelitvijo napak (defektov) v snovi in Weibullovo porazdelitvijo, pa še to za preproste preskusne geometrije (homogeni vzorci, homogena napetost v vzorcu) (9). V tem prispevku je opisana metoda simulacije Monte Carlo porazdelitve robnih razpok v vzorcu in numerične simulacije enoosnega upogibnega preskusa. Na osnovi mnogokratne ponovitve takšne simulacije (vsaka simulacija za en preskusni vzorec) se da napovedati statistično porazdelitev upogibne trdnosti vzorcev pri dani preskusni geometriji. Pri tem lahko odgovorimo na dvome, ali imamo Weibullovo porazdelitev merjene mehanske veličine tudi za kompleksnejše geometrije, na primer pri večplastnih kompozitih. 2 TEORETIČNA OSNOVA ZLOMA KRHKEGA MATERIALA Resnična trdnost keramike je več velikostnih redov manjša od teoretične trdnosti, izpeljane iz medatom-skih sil pri idealnem materialu brez napak. Vzrok za to so različne vrste napak (defektov): (mikro)razpoke, pore, fazne meje itd. Zaradi geometrijskih razlogov je lahko mehanska napetost ob defektu veliko večja kot drugod v snovi, torej večja tudi od nominalne napetosti zaradi zunanjih sil na vzorec. Povečanje napetosti lahko izračunamo z linearno elastično teorijo deformacije za defekte preprostih geometrijskih oblik, npr. ozke razpoke in okrogle pore. Pogoj za katastrofalno širjenje razpoke, ki se je začela ob nekem defektu, lahko izpeljemo z minimizacijo proste energije ali pa neposredno s primerjavo dejanske napetosti ob konici razpoke in kritične napetosti. Obakrat pridemo do neenačbe: K > Kr (1) kjer je K faktor intenzitete napetosti ("stress intensity factor" v literaturi, na kratko imenovanega faktor K), Kc pa kritični faktor intenzitete napetosti, ali drugače, zlomna žilavost. Zlomna žilavost Kc je inherentna lastnost snovi, povezana z medatomskimi vezmi in jo lahko izrazimo z Youngovim modulom in površinsko energijo pri prelomu. Faktor intenzitete pa je odvisen od napetosti (daleč od defekta oz. lokalne napetosti v materialu, če bi bil brez defektov, npr. vsiljene napetosti pri upogibnem preskusu) in dolžine razpoke a, glede na katero ta faktor izračunamo: K = Yo4a (2) kjer je Y brezdimenzijski geometrijski faktor, odvisen od oblike in orientacije defekta (razpoke). Za katastrofalno širjenje razpoke mora biti izpolnjena neenačba (1), kar pomeni, da lahko iz velikosti največjih razpok (defektov) izračunamo, pri kolikšni napetosti se bo material zlomil. Čeprav je v materialu veliko defektov, se bo zlom začel pri enem od največjih. Ker je keramika veliko bolj občutljiva za natezne kot za tlačne napetosti, računamo faktor K le tam, kjer je napetost natezna (pozitivna). Če so v materialu tudi preostale napetosti, npr. termične napetosti zaradi različnih temperaturnih razteznostnih koeficientov komponent v kompozitu, jih moramo prišteti k napetostim zaradi zunanjega obremenjevanja vzorcev. Enačba (2) velja, če je napetost v materialu homogena, npr. pri direktnem nateznem preskusu (če izvzamemo njeno koncentracijo ob defektih), sicer pa si pomagamo pri računu faktorja K z integralom produkta nehomogene napetosti in posebne utežne funkcije (ki je odvisna od oblike razpoke) po dolžini razpoke. Podrobnosti najdemo npr. v člankih (10-12). Pri upogibnem preskusu se mehanska napetost linearno spreminja vzdolž dožine razpoke, ki je pravokotna na spodnjo ploskev. Shema 4-točkovnega enoosnega upogibnega preskusa je prikazana na sliki 1. Največja natezna napetost, označimo jo s amax, se pojavi na spodnji ploskvi vzorca, v območju L/2 - LyJ2 < ^ < L/2 + Lin/2, y = 0, to je v razmiku notranjih podpor Lin. Od tega območja se mehanska napetost linearno spreminja v x- in y-smeri, pri plastnatih kompozitih pa ima na mejah plasti preskoke zaradi različnih Youngovih modulov (13). Pri y = t/2 (t je debelina vzorca) in x = 0 ali x = L (lega zunanjih podpor v razmiku L) je ta napetost enaka nič. Pričakujemo, da se bo zlom začel v območju največje mehanske napetost^max, in to napetost, ko so sile ravno dovolj velike za zlom, imenujemo navadno upogibna trdnost materiala in jo označimo Pojavijo se lahko tudi izjeme: Če je neka izrazito velika napaka (defekt) zunaj območja največje napetosti, npr. desno od tega območja, se zlom lahko začne tam, in to je treba upoštevati pri računu dejanske upogibne trdnosti. Poleg tega je treba pri natančnem računu upoštevati tudi to, da je ob konici razpoke (y > 0) mehanska napetost malo manjša kot pri y = 0. V osnovnem približku ne upoštevamo robnih nehomogenosti napetosti in je pri navadnih preskusnih geometrijah ta približek upravičen. Napetost je tudi neodvisna od koordinate z (pravokotno na ravnino slike). Slika 1: Shema 4-točkovnega upogibnega preskusa. Če je Lin = 0, je to 3-točkovni upogibni preskus. Prikazan je 3-plastni simetrični kompozit z zaporedjem izmenjajočih se plasti A/B/A. Na sliki 1 so ponazorjene tudi ozke robne razpoke, naključno porazdeljene na spodnji ploskvi (kratke črtice), naključne pa so tudi njihove velikosti. V tem prispevku je opisana njihova simulacija Monte Carlo. Realni vzorci imajo še druge vrste defektov, vendar nam omejitev opisa na robne razpoke daje vse bistvene poteze porazdelitve zlomnih sil/napetosti. 3 WEIBULLOVA PORAZDELITEV Dvoparametrično Weibullovo porazdelitev lahko ekvivalentno opišemo z eno ali drugo od dveh porazdelitvenih funkcij: N m-1 f ( o u) = ■ m o o v0 0y exp m o v 0 0 y F (o u) = f( o )do = 1 - exp / \ m \ ^ o u ' v o 0 y (3a) (3b) Spremenljivka o^ = o^ax je upogibna trdnost danega merjenca (nominalna največja napetost na spodnji - natezni strani vzorca pri upogibnem preskusu), ki jo lahko izračunamo iz geometrijskih parametrov in zlomne sile. Ta trdnost se razlikuje od vzorca do vzorca, zato govorimo o porazdelitvi upogibne trdnosti. Funkcija f je verjetnostna porazdelit-vena funkcija z običajnim pomenom: f(ou)dou je verjetnost, da bo ležala izmerjena upogibna trdnost v majhnem intervalu ^u^u+^u). Funkcija F je kumulativna porazdelitvena funkcija, ali drugače, zlomna verjetnostna funkcija. Njen pomen je naslednji: F^„) pomeni verjetnost, da bo izmerjena upogibna trdnost manjša od vrednosti a„. Medtem ko ima fukcija f maksimum, funkcija F monotono narašča od 0 do 1, ko raste njen parameter od 0 do neskončne vrednosti. Seveda, verjetnost, da se bo material prej ali slej zlomil, če napetost poljubno povečujemo, je 1. Weibullova parametra sta Weibullov modul m in karakteristična trdnost 0. Z njima je mogoče izraču- o nati srednjo vrednost upogibne trdnosti in njeno standardno deviacijo: < ou >= ou r(1 + 1 / m) (4a) SD( o u) = o (4b) kjer j^ gama funkcija. Dobra keramika ima visok Weibullov modul, m > 10. Pri tolikšnih modulih je pričakovana vrednost trdnosti ^„> odvisna predvsem od karakteristične trdnostm0 (^u> je malo manjša kot a0) in zelo malo od m, saj j^(1+1/m^ 1. Standardna deviacija pa je močno odvisna od obeh parametrov; čim večji je m, tem manjša je SD^„). Visok Weibullov modul tako daje keramičnemu izdelku večjo zanesljivost, saj če je porazdelitvena funkcija trdnosti ožja, je s tem tudi veliko manjša verjetnost, da bo izdelek odpovedal pri majhnih obremenitvah. Že za relativno skromne vrednosti Weibullovega modula m je Weibullova porazdelitvena funkcija navidezno takšna kot Gaussova funkcija z enako srednjo vrednostjo in standardno deviacijo. Vendar pa se Weibullova porazdelitev v vsaj dveh stvareh bistveno razlikuje od Gaussove. Prvič, je nesimetrična glede na svojo srednjo vrednost in njeni "repi" (tam, kjer se naključna spremenljivka bistveno razlikuje od svoje srednje vrednosti), čeprav navidez nepomembni, se razlikujejo od Gaussovih "repov". To pomeni, da je verjetnost, da se izdelek zlomi pri relativno majhnih obremenitvah, drugačen pri Weibullovi porazdelitvi kot pri Gaussovi. Drugič, Weibullova porazdelitev vodi do znane kvantitativne ocene za vpliv velikosti vzorca na njegovo srednjo upogibno trdnost, čeprav bi na prvi pogled morala biti od geometrije vzorcev in preskusov odvisna le zlomna sila, ne pa iz nje izračunana trdnost. Velja naslednja enačba za odvisnost karakteristične trdnosti od preskušanega volumna vzorca: < o 01 > V 1 \ — V1 (5) vpliva geometrije na izmerjene vrednosti trdnosti keramičnih vzorcev je treba natančno navesti geometrijske podatke pri upogibnem preskusu. Weibullova statistika teoretično izhaja iz določene porazdelitve velikosti napak v materialu. Če imamo samo en tip napak, npr. robne razpoke, je mogoče pokazati '3), da je v primeru naslednje porazdelitve velikosti a napaka: f (a) = /0 a- (6) Pri tem sta V1 in V2 efektivna obremenitvena volumna dveh serij preskušanih vzorcev iz istega materiala, a01 in a02 pa ustrezna karakteristična parametra Weibullove porazdelitve. Enako razmerje kot za karakteristični trdnosti velja tudi za srednji vrednosti trdnosti obeh serij. To pomeni, na primer, če je Weibullov modul keramičnega materiala 10, prva serija vzorcev pa ima pri 3-točkovnem preskusu dvakrat večji razmik med zunanjima podporama (2-krat večji preskusni volumen) kot druga serija, bo povprečna vrednost merjene trdnosti pri prvi seriji manjša kot pri drugi, in sicer za faktor 21/10 = 1,07, to je razlika 7 %. Iz enačbe (5) je razvidno tudi, da so trdnosti vzorcev pri večjem Weibullovem modulu manj občutljive za razlike v preskusnih geometrijah. Zaradi pri čemer je f verjetnostna porazdelitvena funkcija, f0 normalizacijski koeficient, q pa značilni eksponent, porazdelitev trdnosti res Weibullova z Weibullovim modulom m = 2(q - 1). Če porazdelitev velikosti napak ne ustreza funkciji (6), je treba za vsak primer posebej preveriti, ali je porazdelitev trdnosti vsaj približno Weibullova ali ne, bodisi z analitičnim verjetnostnim računom ali pa numerično '9). Dalje, če imamo več kot en tip napak v materialu, ki vse bistveno vplivajo na njegovo trdnost, potem preprosta 2-parametrična Weibullova porazdelitev (3) ne zadošča in je treba za opis porazdelitve trdnosti vzorcev obravnavati multimodalno Weibullovo porazdelitev. Iz relativno majhnega števila izmerjenih/simuliranih vrednosti trdnosti se da z uporabo funkcije (3b) in različnih postopkov prirejanja določenih verjenosti F danim podatkom a„ oceniti vrednosti Weibullovih parametrov m i^0 '2). Zavedati pa se je treba, da je to le ocena in da se lahko predvsem dejanska vrednost modula m precej razlikuje od njegove ocene. Vsekakor velja splošno pravilo statistike: čim večje je število izmerjenih vrednosti naključne spremenljivke, tem zanesljivejša je ocena parametrov porazdelitve. Standardi priporočajo vsaj 30 meritev za zadovoljivo zanesljivost ocene Weibullovih parametrov (14). Ena od metod ocenjevanja Weibullovih parametrov je linearna regresija (LR) (2). 4 SIMULACIJA MONTE CARLO PORAZDELITVE RAZPOK IN WEIBULLOVE PORAZDELITVE Z uporabo računalniškega generatorja naključnih števil med 0 in 1 in ustrezno funkcijsko transformacijo lahko generiramo kakršno koli porazdelitveno funkcijo, npr. porazdelitev velikosti razpok (6). Numerični postopek, ki smo ga sestavili za vedenje vzorcev z robnimi razpokami pri upogibnih preskusih, gre takole. Za vsak namišljeni vzorec posebej se pri upogibnih preskusih z generatorjem naključnih števil generira določeno število razpok z naključnimi koordinatami x in velikostmi a (slika 1). Potem se postopoma povečuje sila F pri preskusu in za vsako razpoko posebej preveri, ali ni morda že dosežen pogoj (1) za začetek njenega širjenja. Če se neka m razpoka pri dani sili začne širiti, je treba pogoj njenega širjenja po majhnih prirastkih njene dolžine spremljati do preloma, saj se lahko pri zapletenih geometrijah, npr. pri plastnatih kompozitih, zgodi, da se razpoka razširi samo do določene plasti, potem pa se njeno širjenje ustavi. Ko se pri dani sili vsaj ena od razpok razširi do končnega zloma vzorca, se ugotovi ustrezna mehanska napetost na spodnji strani vzorca in se opredeli kot upogibno trdnost za dani vzorec. Ta postopek je treba ponoviti za celo serijo vzorcev (N vzorcev, npr. N = 100), da lahko izračunamo parametre porazdelitve trdnosti, da ugotovimo, ali je porazdelitev sploh Weibullova. To nam pove korelacijski koeficient pri ocenjevalni metodi LR: Če je ta koeficient večji kot 0,9, lahko trdimo, da metoda Monte Carlo metoda dobro simulira Weibullovo porazdelitev. Zaradi statističnih deviacij bomo dobili pri drugi seriji 100 vzorcev nekoliko drugačen rezultat za ocenjene Weibullove parametre kot pri naslednjih serijah. Zato ponovimo mnogokrat simulacijo na različnih serijah (seveda pod enakimi pogoji), Np-krat. Če vzamemo npr. NP = N = 100, pomeni to meritve na 100 serijah po 100 vzorcev, to je 10 000 simuliranih upogibnih preskusov. Ker je NP = 100, dobimo 100 nekoliko različnih vrednosti Weibullovih parametrov. Smiselno je povprečiti 100 vrednosti za končno oceno parametrov, zanimivo pa je omeniti, da je rezultat bolj pravilen, če povprečimo njihove logaritme. S tem da smo dobili 100 različnih vrednosti parametrov, smo lahko ugotovili tudi zanesljivost njihove končne ocene. 5 REZULTATI Metodo smo najprej preskusili na homogenem vzorcu iz aluminatne (Al2O3) keramike. Iz literature(12) smo vzeli podatek za njeno zlomno žilavost KC = 3,8 MPa m1'2. Zaradi priročnosti smo razmik med zunanjima podporama normalizirali na L = 1. Porazdelitev upogibnih trdnosti se v vsakem primeru odlično ujema z Weibullovo porazdelitvijo in je močno odvisna od eksponenta q v enačbi (6). Ne le, da parameter m v večini primerov ustreza enačbi m = 2(q — 1), tudi karakteristični parameter 0 je močno odvisen od q. 3-točkovni preskus (Lin = 0) daje enak Weibullov modul, a večje vrednosti karakterističnega parametra kot 4-točkovni, kar je v skladu z vplivom efektivnega volumna vzorca. To je zato, ker je povprečna mehanska napetost pri isti sili F po volumnu vzorca manjša pri 3-točkovnem preskusu. Dvakrat večji preskusni volumen pri sicer enakih pogojih smo simulirali z 2-krat večjim številom simuliranih robnih razpok in ugotovili odlično ujemanje z enačbo (5). Ugotovili smo tudi, da se pri primernem plastnatem kompozitu, kot je A/AZ/..../A, kjer si izmenično sledijo plasti iz Al2O3 (oznaka A) in plasti kompozita Al2O3-ZrO2 (aluminijev in cirkonijev oksid, oznaka AZ), lahko bistveno povečata oba Weibullova parametra. To pomeni večjo efektivno trdnost kompozita v primerjavi s homogenim vzorcem Al2O3. Vzrok za povečanje trdnosti so tlačne preostale termične napetosti v plasteh A (13,15). Čeprav je bil originalno v porazdelitvi velikosti razpok (6) mišljen interval njihovih velikosti od 0 do 00, smo vzeli v svojih simulacijah omejitev amn < a < amax in seveda funkcijo (6) ustrezno normalizirali. Za to smo imeli tri razloge: 1) neskončni interval začetnih dolžin razpok fizikalno ni smiseln, 2) zaradi nume-ričnih razlogov mora biti interval končen, 3) z intervalom (amin, a^^i) se da lepo prilagoditi karakteristični parameter 0 eksperimentalnim vrednostim, ne da bi bistveno spremenili Weibullov modul. Vendar pa lahko prevelika vrednost a^in precej spremeni Weibullov modul m glede na njegovo pričakovano vrednost 2(q - 1), zato je treba iskati v parametrizaciji kompromise, tudi glede na časovno zahtevnost simulacij. Nekaj rezultatov za Al2O3 prikazujeta tabeli 1 in 2. Tabela 2 vsebuje realnejše vrednosti karakterističnega parametra kot tabela 1, vendar pa se pri q = 3,5 precej razlikuje od pričakovane vrednosti 5. Med drugim je zanimiva primerjava med 4-točkovnim (4T, razmik med notranjima podporama je polovica razmika med zunanjima podporama) in 3-točkovnim (3T) preskusom. Pri še večjem razmiku med notranjima podporama (ni prikazano v tabelah), npr. Lin = 0,8 L, se karakteristična trdnost v skladu s pričakovanji še zmanjša. Tabela 1: Izračunana Weibullova parametra m and Oo iz simulacij Monte Carlo za Al2O3 pri upogibnih preskusih. Vrednost (m) v oklepaju v prvem stolpcu ustreza enačbi m = 2(q - 1). Drugi parametri simulacije: amin = 0,5 ^m, amax = 50 ^m, 5000 razpok, debelina vzorca 100 ^m, N = NP = 100. q(m) Lin/L mOc (msp-mZG)