Računalniško podprta identifikacija temperaturne odvisnosti snovnih lastnosti Computer Aided Identification of the Temperature Dependence of Material Properties B. Štok1, P. Koc, FS Ljubljana Prejem rokopisa - received: 1995-10-04; sprejem za objavo - accepted for publication: 1995-12-22 Današnje stanje razvoja matematičnega in numeričnega modeliranja omogoča obravnavo večine tehničnih problemov na osnovi numerične analize in računalniške simulacije. Verodostojnost takšnih simulacij je pogojena v veliki meri s stopnjo poznavanja snovnih lastnosti, ki obravnavani problem karakterizirajo. Še posebej je pomembna temperaturna odvisnost, pa naj gre za toplotne, reološke ali elektromagnetne lastnosti snovi. V prispevku obravnavamo računalniško podprto identifikacijo temperaturne odvisnosti snovnih lastnosti, zasnovano na metodologiji reševanja inverznih problemov. Ključne besede: matematično modeliranje, inverzni problem, snovna identifikacija The actual state of the art in mathematical and numerical modelling has enabled most of technical problems to be treated by means of numerical analysis and computer simulation. The reliability of such simulations is affected considerably by the degree to which the material properties characterizing the considered problem are specified. In technical applications it is the temperature dependence of thermal, rheologic and electromagnetic material properties which is of particular interest. The computer aided identification of the temperature dependent material properties, vvhich is based on the inverse problem solution methodology, is presented in the paper. Key words: mathematical modelling, inverse problem, material characterization 1 Uvod Ugotavljanje velikosti snovnih lastnosti temelji izključno na matematičnem vrednotenju ter obdelavi eksperimentalno pridobljenih odzivov, pri čemer mora biti za nedvoumno ter enolično identifikacijo raziskovane fizikalne lastnosti eksperiment ustrezno načrtovan ter izveden. Zal pa kompleksnost odvisnosti, ki naj bi jo iz eksperimenta identificirali, pogojuje zahtevnost in število potrebnih eksperimentov, s tem pa tudi ceno izvedbe eksperimenta. Da bi enolično identificirali neko odvisnost, temelji večina tako predpisanih eksperimentov na zagotavljanju konstantnosti vplivnih veličin med eksperimentom. Takšen način pa zahteva pri ugotavljanju kompleksnejših odvisnosti izvedbo množice preizkusov v različnih razmerah. V prispevku je obravnavana računalniško podprta identifikacija temperaturne odvisnosti snovnih lastnosti, ki temelji na računalniški simulaciji laboratorijskega eksperimenta ter rešitvi inverznega robnega problema. S tem je podano izhodišče za popolnejšo in hitrejšo karakterizacijo snovi, z možnostjo istočasne sprostitve opisanih omejitev pri izvedbi eksperimenta. ' Prof.dr. Buriš ŠTOK Univerza v Ljubljani - Fakulteta 7.a strojništvo Laboratorij za numcrično modeliranje in simulacijo v mehaniki 10(X) Ljubljana. Aškerčeva 6 2 Identifikacija temperaturne odvisnosti toplotnih lastnosti Pri obravnavanju problemov, katerih časovni razvoj sledi spreminjajočemu se temperaturnemu polju, je poznavanje temperaturne odvisnosti snovnih lastnosti odločujoče za verodostojno analizo. V nadaljevanju se bomo omejili na identifikacijo snovnih lastnosti, ki vplivajo na prevod toplote v trdninah. 2.1 Enačba prevoda toplote Spreminjanje temperaturnega polja T = T(xbt) v trdni snovi določa diferencialna enačba prevoda toplote 3T p c — = div(k grad T) + Q (1) ot pri čemer je prostorska in t časovna koordinata. V enačbi so nadalje T temperatura, p specifična gostota, c specifična toplotna kapaciteta, k toplotna prevodnost in Q polje prostorsko porazdeljenih toplotnih izvirov. V smislu rešitve direktnega robnega problema, s predpisanim časovnim spreminjanjem robnih pogojev ter prostorske porazdelitve toplotnih izvirov, daje diferencialna enačba ob poznanem gradivu s snovnimi lastnosti p, c in k, ki so v splošnem temperaturno odvisne, ter poznanih začetnih pogojih časovno spreminjanje temperaturnega polja T = T(xk,t). V smislu identifikacije v enačbi (1) prisotnih snovnih veličin lahko glede na zgradbo enačbe, v kateri se p in c pojavljata v obliki produkta, k pa v povezavi z diferencialnim operatorjem, ugotovimo, da je v postopku identifikacije mogoče enolično določiti le odvisnosti toplotne prevodnosti k in toplotne kapacitete p c. Enolična identifikacija specifične gostote p in specifične toplotne kapacitete c hkrati ni možna. Zato bomo v nadaljevanju privzeli, da je temperaturna odvisnost specifične gostote p = p(T) poznana, iskali pa bomo temperaturni odvisnosti specifične toplotne kapacitete c = c(T) in toplotne prevodnosti k = k(T). 2.2 Opis inverzne metode identifikacije Računalniško podprta identifikacija temperaturne odvisnosti analizirane snovne lastnosti je zasnovana na privzetju funkcijsko opredeljene odvisnosti, s temperaturo kot neodvisno spremenljivko. Ne glede na izbiro funkcijskega predpisa je v postopku računske identifikacije potrebno določiti velikost koeficientov, ki funkcijsko odvisnost opredeljujejo, skladno z naravno odvisnostjo. Zaradi tega lahko funkcijske koeficiente c i, C2,..., c„ pri opisu funkcijske odvisnosti c = c(T) ter koeficiente kj, k2,..., km pri opisu funkcijske odvisnosti k = k(T) opredelimo kot snovne konstante. Snovne konstante C/, C2,..., c„ in ki, k2,..., km so v obravnavanih funkcijskih odvisnostih z vidika identifikacije neodvisne spremenljivke, katerih velikost naj bi v procesu identifikacije določili. Združimo jih v vektor neznank pT = [ci,...,km], pri čemer označuje T transponiranost vektorja p. Ker gre za inverzni problem, katerega reševanje je enako načinu reševanja ek-stremalnih problemov2, bomo v procesu računske identifikacije iskali optimalno vrednost vektorja neznank p z ozirom na mi-nimizacijo odstopanja med računskim in eksperimentalnim odzivom. Ciljno funkcijo, katere minimum iščemo, definirajmo kot vsoto kvadratov v matrični obliki S(p) = [Y - U(p)]TW[Y - U(p)] + a((Hp)THp) (2) kjer je Y vektor na opazovanem sistemu izmerjenih temperatur in V (p) vektor z numeričnim modelom izračunanih temperatur. Slednjega, ki je implicitna funkcija parametrov sistema p, izračunamo na osnovi analize direktnega problema. Napake meritev vključuje na statistični način utežna matrika medtem ko je vloga regularizacijskega člena z utežnim koeficientom (X, da zmanjšuje prevelike spremembe vrednosti parametrov sistema, ki so navadno značilne za začetne korake identifikacije. Z regularizacijsko matriko H izvajamo nad parametri sistema regularizacijo 0., 1. ali 2. reda, pri čemer določa utežni koeficient a delež regularizacijskega člena v iznosu ciljne funkcije. Potrebni pogoj za ekstrem ciljne funkcije S(p), s katerim izrazimo zahtevo po ničnosti njenih prvih odvodov po parametrih p, zapišemo z enačbo VpS(p) = -2XT(p)W[Y - U(p)] + 2aHTHp = 0 (3) Pri tem označuje parcialni odvod po komponentah vektorja p in X(p) občutljivostno matriko XT(p) = VpUT(p) (4) Skladno z Gauss-Newtonovo metodo minimizacije ciljne funkcije (2) izvedemo nadalje linearno aproksimacijo vektorja U(p) v (v+1)-vi iteraciji U(pv+1) = Uv+1 = IT + X(pv)[pv+1 - pv] (5) pri čemer pomeni v števec iteracij. Upoštevajoč dano aproksimacijo, sledi iz enačbe (4) sistem enačb, iz katerega ob poznanem vektorju pv ter njemu pripadajočem odzivu U(pv) izračunamo naslednji približek pv+l (XT(pv)WX(pv) + aHTH)[pv+1] = = XT(pv)W(X(pv)[pv] + [Y - U(pv)]) (6) Zapis enačbe (6) predpisuje iteracijski način reševanja inverznega problema. V vsakem iteracijskem koraku moramo izvesti analizo direktnega problema z izračunom vektorja U(pv), ob tem pa še analizo občutljivosti tega vektorja glede na spremembo vektorja pv. Rezultat iteracije Gauss-Newtonove metode je nov vektor parametrov sistema pv+1, ki vodi k minimi-zaciji obravnavane ciljne funkcije. 3 Numerična implementacija modela Izhodišče za računalniško podprto identifikacijo predstavlja množica med potekom realnega eksperimenta izmerjenih vrednosti. Da bi bil posneti odziv čim bolj verodostojen, s tem pa tudi računska identifikacija sama, je potrebno v eksperimentu zagotoviti ustrezno krajevno in časovno zajemanje odziva. Ker nam gre v tem prispevku izključno za prikaz metodologije računske identifikacije, bomo realni eksperiment zamenjali z numeričnim. Neznane odvisnosti snovnih lastnosti v realnem eksperimentu bomo v numeričnem eksperimentu nadomestili s povsem določenimi odvisnostmi, nakar bomo pri predpisanih robnih pogojih z rešitvijo direktnega problema poiskali odziv v eksperimentu obravnavanega sistema. V našem primeru je to časovno spremenljivo temperaturno polje, iz katerega bomo za potrebe računske identifikacije izločili končno množico vrednosti, diskretiziranih tako krajevno kot tudi časovno. S takšno diskretizacijo zagotovimo analogijo z realnim eksperimentom, v katerem je merjenje temperature omejeno krajevno na točke z nameščenimi senzorji, časovno pa na ustrezno časovno zaporedje. Problem obravnavane identifikacije, ki ga matematično opredeljujejo enačbe predhodnega razdelka, rešujemo z algoritmom, opisanim v drugem razdelku. Podrobnejši opis posameznih točk tega algoritma je naslednji: a) Definiramo razmeroma preprost realni eksperiment z nezahtevnim vzdrževanjem robnih pogojev ter ga nadomestimo z ekvivalentnim numeričnim eksperimentom. Analizo slednjega izvedemo z metodo končnih elementov, nakar iz dobljene rešitve izločimo nabor odzivov, ki bo rabil za nadaljnjo identifikacijo. V nasprotju z rezultati numeričnega eksperimenta je pri vsakem realnem eksperimentu v izmerjenih vrednostih vselej prisotna napaka meritve. Ob sistemski napaki meritve, ki je navadno znana in jo lahko ustrezno kompenziramo, preostaja v rezultatih meritev še slučajna napaka neznane velikosti. Maksimalna velikost le-te je odvisna od natančnosti senzorjev in merilne verige in je statistično določljiva. Zaradi navedenega je smotrno, s ciljem približanja razmeram pri identifikaciji na osnovi realnega eksperimenta, numerično določeni nabor odzivov ustrezno perturbirati in tako spremenjeni odziv uporabiti v nadaljnjem postopku identifikacije. V nadaljevanju se sklicujemo na tako določeni nabor kot na odziv realnega eksperimenta. Veličina, ki jo nabor vsebuje, je temperatura v točkah, kjer bi bili pri realnem eksperimentu postavljeni temperaturni senzorji. b) Izberemo obliko funkcijskih odvisnosti, s katerimi bomo aproksimirali iskane temperaturne odvisnosti, ter začetni približek za vektor parametrov sistema p, t.j. snovnih konstant ci,..., c„ in ki,..., k,„. Ker je algoritem identifikacije dovolj robusten, problem prevoda toplote pa dovolj dobro pogojen, izbira začetnega približka ne predstavlja problema. c) Računalniško simulacijo realnega eksperimenta izvedemo na numeričnem modelu, pri čemer upoštevamo v dani fazi računske identifikacije veljavni približek temperaturnih odvisnosti. Analizo tako definiranega direktnega problema izvedemo z metodo končnih elementov. Da bi zagotovili podobnost z identifikacijo pri obravnavi dejanskega realnega eksperimenta, sta pri tem krajevna in časovna diskretizacija problema v splošnem različni od diskretizacije, uporabljene v analizi ekvivalentnega numeričnega eksperimenta. Ker je rezultat računalniške simulacije odziv numeričnega modela v opisanih razmerah, je smiselno, da diskretizacija vključuje diskretizacij-ske točke, ki opredeljujejo nabor odzivov realnega eksperimenta. S tem je namreč enolično določen ekvivalentni nabor odzivov simuliranega numeričnega eksperimenta, ki nam rabi v nadaljnjem postopku za presojo kvalitete dosežene stopnje materialne karakterizacije. č) Presojo kvalitete v dani fazi računske identifikacije dosežene stopnje materialne karakterizacije izvedemo posredno s primerjavo razhajanj numeričnega in realnega odziva. Velikost tega razhajanja spremljamo s ciljno funkcijo S(p), katere zgradbo smo definirali z enačbo (2). Teoretično je v primeru, ko sta odziva realnega sistema in numeričnega modela enaka, vrednost ciljne funkcije S(p) = 0, kar predstavlja absolutni minimum dane funkcije. V procesu računske identifikacije je ciljna funkcija vselej pozitivna, pri čemer je njeno odstopanje od nične vrednosti merilo za nadaljevanje oz. konec računskega postopka. d) Merila kvalitete dosežene stopnje materialne karakterizacije v iteracijskem postopku ne moremo postaviti na osnovi spremljanja absolutne napake ciljne funkcije S(p), marveč kvečjemu na osnovi spremljanja relativne napake obravnavane funkcije skozi iteracijski proces. Glede na velikost relativne napake v danem iteracijskem koraku se računski postopek na tem mestu ustavi ali pa nadaljuje. V primeru izpolnjevanja vnaprej predpisane natančnosti, pri čemer z velikostjo le-te neposredno vplivamo na kvaliteto identificiranih odvisnosti, računski postopek identifikacije iskanih snovnih odvisnosti končamo. e) Pri neizpolnjevanju predpisane natančnosti se računski postopek identifikacije iskanih snovnih odvisnosti nadaljuje z ustreznim popravkom snovnih konstant. Novi približek parametrov sistema določimo z Gauss-Newtonovo metodo, skladno 1200 1100 1000 « 900 bc \ 800 700 600 500 400 -•- C —A— p \ - 100 8000 90 7900 BO 7800 ,—, 70 « - 7700 ra 6 b 60ž " 7600 M 50^ - 7500 o. 40 7400 30 7300 20 7200 0 200 400 600 800 1000 1200 T [°C] Slika 1: Toplotne lastnosti preizkušanca Figure 1: Thermal properties of the sample z enačbami (3-6), nakar izvedemo vrednotenje dane spremembe s ponovitvijo računskega postopka od točke c) dalje. 4 Numerični zgled V ponazoritev prikazane metodologije računalniško podprte snovne identifikacije smo izvedli identifikacijo temperaturno odvisnih toplotnih lastnosti mehkega železa5. Temperaturne odvisnosti specifične toplotne kapacitete, toplotne prevodnosti in specifične gostote danega materiala so prikazane na sliki 1. Zaradi razlogov, omenjenih v razdelku 2.1, sta predmet identifikacije samo specifična toplotna kapaciteta c in toplotna prevodnost k. Sklop za izvedbo laboratorijskega eksperimenta, katerega geometrija je prikazana na sliki 2, je sestavljen iz dveh enakih plošč neznanega materiala (mehko železo) ter ploščatega grelnika z močjo 500 W, ki je vstavljen med plošči. Tesno prilegajoča toplotna izolacija (šamot z lastnostmi: c = 840 + 0,322 (T-100) v J/kgK ter k = 0,8 + 0,00022 (T-100) v W/mK. pri čemer je T v °C) preprečuje konvektivni prenos toplote v okolico. Senzorji za merjenje temperature so postavljeni med grelnikom in ploščo (oznaka A na sliki 2) ter med ploščo in izolacijo (oznaka B). Če imamo na voljo več senzorjev, jih razporedimo po ploskvah obeh plošč in s tem povečamo zanesljivost odčitanih meritev. Iz začetnega stanja, ko ima cel sklop grelnik-plošči-izolacija temperaturo okolice T„, segrevamo plošči, vse dokler senzorji na površini izolacije (oznaka C) ne zaznajo dviga temperature. Ker je od takrat dalje možen prenos toplote s konvekcijo, s čemer se robni pogoji v smislu njihovega nadzora bistveno spremenijo, segrevanje ustavimo in prekinemo odvzem podatkov s senzorjev A in B. Opisani laboratorijski eksperiment nadomestimo v našem primeru z analizo ekvivalentnega numeričnega eksperimenta. Da bi zagotovili podobnost z realnim odzivom, perturbiramo numerični odziv skladno z napako meritve, ki izvira iz netočnosti senzorjev. Glede na temperaturno območje ter vrsto senzorjev je določena največja verjetna napaka, ki je Terr = ±max(2,50C; 0.0075AT) , (7) pri čemer \T pomeni temperaturno razliko med toplim in hladnim stikom termoelementa. Model funkcijske odvisnosti, ki smo ga v računski identifikaciji privzeli za aproksimacijo iskanih temperaturnih odvisnosti, je odsekoma linearna funkcija. Pri izbranih dimenzijah Slika 2: Eksperimentalni sklop Figure 2: Experimental arrangement obravnavanega sklopa in dani kvaliteti toplotne izolacije se izkaže, da z enim samim eksperimentom, izvedenim tako, da vključuje celotni temperaturni interval od 0°C do 1300°C, ni mogoče izvesti zadovoljive identifikacije. Zato celotni eksperiment razdelimo na več eksperimentov, razpotegnjenih preko več podintervalov. V tabeli 1 je za izbrana temperaturna območja prikazano število za identifikacijo potrebnih iteracij in končna velikost ciljne funkcije. Začetna približka privzemata konstantno vrednost snovnih lastnosti na celotnem temperaturnem območju od 0°C do 1300°C, za specifično toplotno kapaciteto co = 400 J/kgK ter toplotno prevodnost k0 = 50 W/mK. Zaradi izjemne skokovitosti krivulj v območju od 800°C do 1000°C je ta del diagrama zajet z dodatnim eksperimentom (številka 5). Tabela 1: Rezultati identifikacije Table 1: Results of the identification eksperiment temp.območje (°C) št. iteracij ciljna funkcija 1 0-425 3 126,8 2 400 - 775 5 112,3 3 700 - 1000 4 134,8 4 850 - 1300 5 139,3 5 800 - 980 5 157,2 Rezultati identifikacije so za vseh pet preizkusov združeni in prikazani na sliki 3. S polnimi simboli sta označeni točni krivulji c = c(T) in k = k(T), z votlimi pa krivulje, ki so rezultat identifikacije. V večjem delu diagrama je ujemanje pravih in ocenjenih vrednosti zadovoljivo, na mestih skokovitih prehodov pa je po pričakovanju ujemanje slabše. 5 Sklepne ugotovitve Iz prikazanega numeričnega primera, ki je bil izbran zaradi zapletene temperaturne odvisnosti snovnih lastnosti, je razvidna robustnost metode identifikacije. Velike nezveznosti, ki se kažejo v skoraj stopničasti funkciji v okolici 920°C, pa so 1200 1100 1000 i 900 800 700 600 500 400 ! II —•— C Va h ■ —B— k-oce na "i t \ \ J / V % S* i" 100 90 80 70 £ 60 \ 1-1 50 40 30 20 200 400 600 800 1000 1200 1400 T [°C] Slika 3: Identificirane lastnosti Figure 3: Identified properties velik problem za vsak algoritem. Obravnavana metodologija identifikacije ni uporabna zgolj za snovno identifikacijo, marveč jo je mogoče smiselno uporabiti za pripravo eksperimenta in pravilno vodenje le-tega. Iz analize občutljivosti (razdelek 2.2) je mogoče iz velikosti koeficientov občutljivostne matrike razbrati optimalno lego senzorjev. 6 Literatura 'J. S. Arora: Introduction to Optimum Design. McGravv-Hill Book Company, 1989 2 J. V. Beck, K. J. Arnold: Parameter Estimation in Engineering and Science, John Wiley &. Sons, 1977 3 J. V. Beck, B. Blackwell, C. R. Clair: Inverse Heat Conduction - 111-posed Problems, John Wiley & Sons, 1985 4W. J. Minkowycz, E. M. Sparrow, G. E. Schneider, R. H. Pletcher: Handbook of Numerical Heat Transfer, John Wiley & Sons, 1988 5 F. Richter: Die wichtigsten physikalischen Eigenschaften von 52 Eis-enwerkstoffen, Stahleisen - Sonderberichte, Heft 8, Verlag Stahleisen M.B.H., Diisseldorf, 1973