R. KUNC ET AL.: DOLO^ANJE SNOVNIH LASTNOSTI MATERIALA ZA RA^UNANJE NOSILNOSTI ...
DOLO^ANJE SNOVNIH LASTNOSTI MATERIALA ZA RA^UNANJE NOSILNOSTI VELIKIH KOTALNIH
LE@AJEV Z MEHANIKO PO[KODBE
DETERMINATION OF MATERIAL PROPERTIES FOR THE CALCULATION OF LARGE ROLLING BEARING CARRYING
CAPACITY USING DAMAGE MECHANICS
Robert Kunc1, Ivan Prebil1, Jo`e Korelc2, Toma` Rodi~3, Matja` Torkar4
1Univerza v Ljubljani, Fakulteta za strojni{tvo, A{ker~eva cesta 6, SI-1000 Ljubljana, Slovenija 2Univerza v Ljubljani, Fakulteta za gradbeni{tvo in geodezijo, Jamova cesta 2, SI-1001 Ljubljana, Slovenija
3Univerza v Ljubljani, Naravoslovnotehni{ka fakulteta, A{ker~eva cesta 12, SI-1001 Ljubljana, Slovenija 4In{titut za kovinske materiale in tehnologije, Lepi pot 11, SI-1001 Ljubljana, Slovenija
robert.kunc@fs.uni-lj.si
Prejem rokopisa - received: 2000-12-19; sprejem za objavo - accepted for publication: 2001-03-30
V prispevku je prikazan Pedersenov model mehanike po{kodbe kontinuuma ter eksperimentalno numeri~na dolo~itev snovnih lastnosti materiala 42CrMo4 (ISO 683/1) za razli~na stanja obro~a le`aja velikih dimenzij. Model omogo~a simulacijo elasti~nih in plasti~nih specifi~nih deformacij ter napetosti v odvisnosti od {tevila nihajev obremenitve. Model upo{teva izotropno in kinemati~no utrjevanje oziroma meh~anje ter rast po{kodb materiala z ve~anjem {tevila nihajev. Po{kodbeni model je integriran v numeri~ni model kon~nih elementov, s katerim je narejena simulacija odziva materiala in njena primerjava z eksperimentom.
Klju~ne besede: aksialni le`aj velikih dimenzij, kotalni stik, malocikli~na nosilnost, po{kodba, plastifikacija
The article considers Pedersen's damage model and the experimental-numerical determination of the material parameters for various states of bearing ring material 42CrMo4 (ISO 683/1). The model allows an evaluation of elastic and plastic strain and stress distributions as a function of the number of loading cycles. The model includes isotropic and kinematic hardening or softening and the evolution of damage with the number of cycles. The damage model is implemented into a finite-element code, which has been assessed with respect to results of an experimental test.
Key words: large axial rolling bearings, rolling contact, low-cycle carrying capacity, damage, plasticity
NOMENKLATURA A nepo{kodovani prerez AD po{kodovani prerez D po{kodba
Dc kriti~na vrednost po{kodbe E elasti~ni modul
Lijkl tenzor elasti~nih konstant R izotropno utrjevanje
Ro mejna vrednost izotropnega meh~anja
R∞ trenutna mejna vrednost izotropnega utrjevanja R∞,s mejna vrednost izotropnega utrjevanja
S odpornostna energija po{kodbe Xij sredi{~e napetostnega prostora
X∞(n) mejna vrednost koeficienta kinemati~nega utrjevanja
b koeficient izotropnega utrjevanja/meh~anja f napetostni potencial
k velikost elasti~nega obmo~ja mn Ohno-Wang-ov snovni parameter
p primerjalna specifi~na plasti~na deformacija pd za~etni prag po{kodbe
Λ velikost neproporcionalnosti
γ(n) stopnja kinemati~nega utrjevanja εij tenzor majhnih specifi~nih deformacij
∆λ plasti~ni multiplikator ν Poissonovo {tevilo σ napetost
σy napetost te~enja
σmaxn najve~ja napetost v n-tem nihaju σmaxs najve~ja napetost pri stabilizaciji
1 UVOD
Pri dolo~anju nosilnosti po~asi sevrte~ih aksialnih kotalnih le`ajev velikih dimenzij s povr{insko utrjeno te~ino s standardizirano dopustno trajno deformacijo1ali priporo~ilom, ki dolo~a dopustno podpovr{insko nape- tost 2, nedobimo zadovoljivih rezultatov. Razlog jev neustreznosti meril, ki ne upo{tevajo v celoti zna~ilnosti kotalnega stika velikih kotalnih le`ajev. Tudi pri dopolnjenem merilu dopustnih podpovr{inskih napetosti, ki upo{teva malocikli~ne lastnosti materiala le`ajnega obro~a3, je prera~unana trajna deformacija le`ajne te~ine od eksperimentalno izmerjenih rezultatov tudi do 100%
druga~na.
Dolo~iti `elimo realno nosilnost kotalnega stika velikih kotalnih le`ajev ob predpostavki, da realne razmere zunanjih obremenitev povzro~ajo neenakomer- nost porazdelitve zunanje obremenitve po obodu kotalnih elementov in s tem trajno deformacijo te~ine le`aja oziroma njeno plasti~no deformacijo. Zaradi tega smo se odlo~ili za prera~un nosilnosti kotalnega stika ob upo{tevanju mehanike po{kodbe, pri katerem sprem- ljamo spreminjanjeelasti~nein plasti~nedeformacijeter napetosti in nastanek po{kodb materiala kot funkcijo velikosti in {tevila nihajev obremenitve4,6,11.
2 MEHANIKA PO[KODBE 2.1 Mehanski princip
Po{kodbo si razlagamo kot nastanek mikropovr{in oziroma mikroraztrganin, ki so posledica ve~anja {tevila dislokacij, interakcij med njimi in poru{itev medatom- skih vezi. Poru{itve in interakcije dislokacij povzro~ajo nastanek povr{in, ki so glede na posamezne koordinatne smeri razli~nih velikosti. Novo nastale povr{ine so torej skupni u~inek mikropo{kodb pri volumskem elementu - podobno kot specifi~na plasti~na deformacija, ki je povpre~je ve~jega {tevila zdrsov. Efektivno oziroma dejansko napetost v po{kodovanem materialu dolo~imo s predpostavko (1), da se imenski prerez zmanj{a za velikost po{kodovane povr{ine prereza5:
D M n x A A ( , , )r =δ Dx ⇒
δ σ= σ
− =
− =
− F
A AD D
FA ADA
1 1 (1)
Koncept po{kodbe materiala je vgrajen v konstitu- tivni model za majhne specifi~ne deformacije, ki jih razdelimo na elasti~ni in plasti~ni del:
εij =εije +εijp , (2) kjer je zveza med napetostmi in deformacijami dolo~ena z:
σij = −(1 D L) ijkl⋅εkl . (3) Plasti~ne deformacije dolo~a zakon te~enja:
& &
ε λ
ijp σ
ij
= ∂f
∂ . (4)
Napetostni potencial 6je funkcija napetostnega ten- zorja in utrjevalnih komponent Xij in R ter po{kodbe materiala D. Za izotermno stanje dT/dt = 0 je reolo{ki model oziroma evolucijska ena~ba napetostnega poten- ciala dolo~ena z6,11:
f =ùeq − + =(R k) 0 , (5) ùeq = 3 û ûij ij û =ùij − ijùk
2
1
, ij 3δ ,
ùij =
− − σij
D Xij
1 (6)
Z odvajanjem ena~be (5) dobimo:
& & & & &
f f f
X X f
RR f DD
ij ij
ij
= ∂ ij
∂ + ∂
∂ + ∂
∂ + ∂
∂
σ σ , (7)
kjer je:
∂
∂ =
− f
ij D
ij
σ eq
3 2 1
û ù
( ) ; ∂
∂f = −
R 1
∂
∂ f = XijD
ij eq
3 2
û
ù ; ∂
∂ =
− f
D D
ij ij eq
3 2 1
û û ù
( ) (8)
2.1.1 Kinemati~ni model utrjevanja
Kinemati~no utrjevanje opisuje tenzor Xij, ki dolo~a sredi{~e napetostnega prostora. Upo{tevali smo evolucij- ske ena~be, ki so jih predlo`ili Armstrong in Frederick10 ter Chaboch7,8,14,18.
Xij Xijn
n
=
∑
= ( ) 13 ; X&eq( )n = 2X Xij( )n ij( )n
3 (9)
&( )
Xijn = 2
3γ( ) ( )1 ε ,
( ) ( )
( )& (
n nnn
ij
pl c eqn
n m
ijn
X D X
X X
n
− −
∞
) ( )γ λn & (10) 2.1.2 Izotropno utrjevanje / meh~anje
Material secikli~no utrjujeali meh~a, ~esepri procesu cikli~nega obremenjevanja povr{ina te~enja ve~a ali manj{a. Pri simetri~nem enoosnem cikli~nem obremenjevanju s kontrolo raztezka (konstantna ampli- tuda raztezka) se amplitude napetosti pove~ujejo ali manj{ajo. Pri konstantni amplitudi napetosti pa se amplitude raztezka zmanj{ujejo ali ve~ajo. Zasi~enost cikli~nega utrjevanja je enaka rasti plasti~nega raztezka.
Cikli~na krivuljaσ-εle`i nad monotono krivuljo σ-εali pod njo6.
Velikost povr{ine te~enja opisujeta skalarja R in k. R je spremenljivka, ki opisuje izotropno cikli~no utrjevanje oziroma meh~anje materiala. Za~etne vrednosti pri cikli~nem utrjevanju so: k =σyin R = 0: Pri cikli~nem meh~anju pa so: k = σy - R0 in R = R0. Evolucijska ena~ba za spremenljivko R je podana v obliki:
& ( ( , ) &
R b R= ∞ Λ q −Rλ (11) kjer je b snovni parameter, ki dolo~a stopnjo izotropnega utrjevanja ali meh~anja. Parameter R∞(Λ,q) pa dolo~a mejo izotropnega cikli~nega utrjevanja ali meh~anja11,13.
Vpliv neproporcionalnosti obremenitve na cikli~no utrjevanje ali meh~anje materiala smo privzeli po Benallal-u15. Skalar Λ je uporabljen za oceno velikosti neproporcionalnosti obremenitve. Spremenljivka q pa dolo~a pomik plasti~nega raztezka12.
2.2 Mehanika po{kodbe kontinuuma
Mehanika po{kodbe kontinuuma obravnava po{kod- be materiala pri cikli~nem obremenjevanju. Rast po{kodb materiala lahko eksperimentalno dolo~amo na ve~ na~inov 4,5,6. Pri plasti~no deformiranem materialu
pomenijo `e majhne spremembe elasti~nega modula slabitev materiala zaradi rasti po{kodbe5:
K= −(1 D E.)⋅ (12) Nepovratno rast po{kodbe popi{e evolucijska ena~ba
4,11,16,17, ki upo{teva proporcionalni vpliv efektivne plasti~ne deformacije na spremembo po{kodbe. Evolu- cijska ena~ba za popis (modeliranje) krhke in `ilave po{kodbejepodana v obliki:
&
( ) ( )
D (
S E
eq kk
= eq
+ + −
⋅ ⋅
σ ν ν σ
σ
2
2 2
3 1 3 1 2
3
2 1− 2 ⋅
D p p
) & α( );
α( ) ,
p ,za p p za p p
d d
= ≥
<
1
0 . (13)
Mejo pD, pri kateri za~ne po{kodba nara{~ati, dolo~a akumulirana plasti~na deformacija p6:
& &
p= D
− λ
1 , pD =Max p( (D=0)). (14) Teroteri~no pretrganje materiala nastopi, ko je velikost razpoke enaka velikosti elementa; to je, ko dose`e po{kodba D vrednost 1! V ve~ini primerov pa nastopi pretrganje preostalega dela trenutno, kar povzro~i padec kriti~ne vrednosti DC:
DC Max D p DC
u
= ( ( ))= ≈ −1 σ
σ . (15)
3 EKSPERIMENTALNO DELO
Preizku{anci so bili izdelani iz malolegiranega jekla za pobolj{anje42 CrMo 4 (ISO 683/1; W.Nr. 1.7225).
Jeklo je bilo toplotno obdelano tako, da so imeli preiz- ku{anci po toplotni obdelavi razli~no trdoto. Analiza jekla je bila izvr{ena s kvantometrom ARL.
Natezno-tla~ni in cikli~ni preizkusi so bili izvr{eni na univerzalnem hidravli~nem preizku{evali{~u Instron 1255, kjer smo podatke pridobili z merilnim sistemom National Instruments. Raztezek je izmerjen neposredno na preizku{ancu z dinami~nim merilnikom raztezka
"Dynamic Strain Gauge Ekstenziometer 2620-604",
obremenitev pa z natezno-tla~nim senzorjem sile Lebow
3.
Obliko preizku{ancev smo prilagodili tako, da so bile koncentracije napetosti na prehodih razmeroma majhne.
Glavna os preizku{ancev je bila (slika 1) pre~no na smer valjanja materiala.
3.1 Obremenitveni kolektiv
Najprej so bili narejeni monotoni natezno-tla~ni preizkusi do poru{itve preizku{anca (slika 2) nato pa vrsta cikli~nih enoosnih preizkusov s kontrolo amplitude deformacije za dolo~anje cikli~nih parametrov materiala (preizkusi pri razli~nih vrednostih amplitud deforma- cije). V praksi zadostuje za dobro simulacijo odzivov materiala `e manj{e {tevilo enoosnih cikli~nih preizku- sov11.
Da se izognemo lokalnem zo`enju premera merilne dol`ine, je treba natezno tla~ni preizkus izvesti pri kon-
Slika 1:Preizku{anec Figure 1:Test specimen
Slika 2:Krivulje natezno-tla~nih preizkusov Figure 2:Static tensile and compressive test
∆ε
Slika 3:Dolo~anje koeficientov kinemati~nega utrjevanja (kontrola spec. deformacije:∆ε= 1,0%,εsred= 0)
Figure 3:Determination of the kinematics hardening coefficients (strain control:∆ε= 1.0%,εmean= 0)
stantni amplitudi raztezka. Mehanske lastnosti materiala so bile izmerjene pri temperaturi 20 °C.
4 REZULTATI
Kemijska sestava uporabljenega jekla je podana v tabeli 1.
4.1 Dolo~anje snovnih parametrov
Obravnavani po{kodbeni model vsebuje 23 snovnih parametrov. Ve~ino parametrov smo dolo~ili neposredno z eksperimentom (tabela 2), drugepa po inverzni metodi, kjer smo parametre spreminjali tako, da je bilo ujemanje med numeri~nimi in eksperimentalnimi rezultati ~im bolj{e.
4.1.1 Dolo~anje monotonih parametrov materiala Osnovne parametre materiala, kot so elasti~ni modul - E, Poisson-ovo {tevilo ν in napetost te~enja σy, smo dolo~ili s standardnimi monotonimi enoosnimi natezno- tla~nimi preizkusi (slika 2).
4.1.2 Dolo~anje malocikli~nih parametrov izotropnega in kinemati~nega utrjevanja
Parametre za opis kinemati~nega utrjevanja smo dolo~ili s cikli~no krivuljo σ-ε, dobljene s kontrolo
amplitudedeformacije. Za razli~neamplitudedefor- macije ∆ε dobimo po stabilizaciji histerezne zanke ustrezneamplitudeplasti~nespec. deformacije∆εp in amplitudenapetosti∆σ(Slika 3).
∆σ
2 =(R∞,s + +k X) ∞⋅tanh p
+
1 1
γ ∆2ε
+ ⋅
+ ⋅
∞ ∞
X2 3 p X3 3 p
2 2
tanh γ ∆ε tanh γ ∆ε (16) Parametri materiala γ1, γ2, γ3, X∞1, X∞2 in X∞3 so dolo~eni s prileganjem krivulje (ena~ba 16) z ekspe- rimentalno dobljenim rezultatom. Isto~asno se dolo~i tudi parametra R∞,s in k. Iz preizkusa je viden trend zmanj{evanja maksimuma in minimuma napetosti v odvisnosti od {tevila nihajev, kar pomeni, da se material meh~a (slika 6). Ohno-Wang faktorjem1, m2 in m3
dolo~imo s korelacijo simuliranega odziva materiala z eksperimentalnim. Pri tem uporabimo nesimetri~en enoosni cikli~ni preizkus z lezenjem oziroma z zmanj{e- vanjem srednje vrednosti napetosti (slika 8).
Stopnjo izotropnega utrjevanja ali meh~anja dolo~a parameter b. Dolo~imo ga s prileganjem vrednosti iz ena~be (17) glede na preizkuse s kontrolo amplitude spec. deformacije∆ε(slika 4):
Tabela 1:Kemijska sestava jekla 42CrMo4 (v %)
Tabele 1:Chemical composition of steel 42 CrMo 4 (values in %)
C Si Mn P S Cr Ni Mo Cu Al Sn
42CrMo4 0,43 0,26 0,65 0,015 0,021 1,07 0,19 0,16 0,16 0,021 0,006 Tabela 2:Parametri materiala po{kodbenega modela za jeklo 42 CrMo 4, ki je bilo pobolj{ano na 462 HV
Tabele 2:Material damage model parameters of steel 42 CrMo 4, tempered to 462 HV
E σy ν b R∞,S R0 γ1 γ2 γ3 X∞1 X∞2 X∞3 m1 m2 m3 S pd Dc
195000 650 0,3 5 0 380 190 2000 15000 900 250 200 5 5 5 120 0 0,2
∆ε
∆ε
∆ε
∆ε
Slika 4:Dolo~anje koeficientov izotropnega meh~anja Figure 4:Determination of the isotropic hardening coefficient
Slika 5:Dolo~anjefunkcijerasti po{kodbe(kontrola spec. defor- macije:∆ε= 1,0%,εsred= 0)
Figure 5: Determination of the damage growth as a function of accumulated plastic strain (strain control:∆ε= 1.0%,εmean= 0)
σ σ
σmax σmax ε
max max
exp( )
N s
b p N
−
− 11 = −1 − ⋅2 ∆ ⋅ . (17) 4.1.3 Dolo~anje po{kodbenih parametrov pri
malocikli~ni obremenitvi
Za opis po{kodbe potrebujemo tri materialne para- metre: pd, S in Dc. Rast po{kodbev odvisnosti od nihanja obremenitve dolo~a sprememba elasti~nega modula materiala. Povpre~no vrednost parametra S dolo~imo pri razli~nih obremenitvenih kolektivih preizkusa in nape- tostihσ ter rasti po{kodbe v odvisnosti od akumulirane plasti~nedeformacijep6(slika 5):
S
E D dD
dp
=
− σ2 2 1( )2
, p
i N
ip
=
∑
=2
1
∆ε . (18)
5 DISKUSIJA REZULTATOV
Po{kodbeni model smo vgradili v ra~unalni{ki program, ki temelji na metodi kon~nih elementov 19,20. Za razvoj programa smo uporabili splo{en simbolni na~in19,20, ki omogo~a izdelavo u~inkovite, kvadrati~no konvergentne, implicitne numeri~ne sheme tudi v primeru velikega {tevila evolucijskih ena~b in povezanih nelinearnih sistemov. S programom smo analizirali eksperimente (tabela 2).
Iz primerjave numeri~nih in eksperimentalnih rezultatov je opazno zelo dobro ujemanje maksimalne in minimalne napetosti v odvisnosti od {tevila nihajev obremenitve (sliki 6 in 8) kakor tudi velikost in obliko histereznih zank (sliki 7 in 9). Zmanj{evanje amplitud in srednje vrednosti napetosti na za~etku obremenjevanja (do stabilizacije histereznih zank) dobro opi{eta simu-
∆ε ε
Slika 7:Potek histerezσ-ε(kontrola spec. deformacije:∆ε= 1,0%, εsred= 0)
Figure 7:Comparison ofσ-εhysteresis (strain control:∆ε= 1.0%, εmean= 0)
∆ε ε
Slika 9:Potek histerezσ-ε(kontrola spec. deformacije:∆ε= 1,4%, εsred= 0,5%)
Figure 9:Comparison ofσ-ε hysteresis (strain control:∆ε= 1.4%, εmean= 0.5%)
∆ε ε
Slika 8: Potek amplitud napetosti (kontrola spec. deformacije:
∆ε= 1,4%,εsred= 0,5%)
Figure 8:Peak stress vs. number of cycles (strain control:∆ε= 1.4%, εmean= 0.5%)
∆ε ε
Slika 6: Potek amplitud napetosti (kontrola spec. deformacije:
∆ε= 1,0%,εsred= 0)
Figure 6:Peak stress vs. number of cycles (strain control:∆ε= 1.0%, εmean= 0)
lacija izotropnega meh~anja in kinemati~nega utrjevanja.
Po stabilizaciji histereznih zank oziroma po zmanj{anju vpliva kinemati~nega in izotropnega utrjevanja ima glavni vpliv na napetostno-deformacijsko stanje rast po{kodbemateriala. Ve~anjepo{kodbepovzro~i rahlo zmanj{evanjeamplitudenapetosti vsedo kon~ne poru{itve.
6 SKLEPI
Za dolo~anjerealnenosilnosti po~asi vrte~ih se aksialnih kotalnih le`ajev velikih dimenzij s povr{insko utrjeno kotalno te~ino potrebujemo prera~un nelinear- nega napetostno-deformacijskega stanja v le`ajnem obro~u. Numeri~ni prera~un z upo{tevanjem mehanike po{kodb zelo dobro opi{e nelinearni odziv materiala na malocikli~no obremenjevanje tako po velikosti amplitud deformacije in napetosti kakor tudi po obliki histereznih zank.
Z mehaniko po{kodbe `elimo izra~unati velikost podpovr{inskih deformacij in napetosti kakor tudi po{kodbe osnovnega in utrjenega materiala na meji povr{insko kaljene plasti. S tak{no analizo nameravamo dolo~ili kriti~na mesta in {tevilo obremenitvenih ciklov do nastanka po{kodbe ter dobiti oceno o hitrosti {irjenja po{kodb.
Za modeliranje realnih razmer kotalnega stika je treba dolo~iti po{kodbene parametre za celoten prerez le`ajnega obro~a. To pomeni dolo~itev parametrov za povr{insko kaljeno plast te~ine in tudi zmanj{anje napetosti v osnovnem materialu obro~a le`aja pod kaljeno te~ino.
7 LITERATURA
1ISO 76, Rolling bearings - Static load ratings, 1987
2Pallini, R. A.; Sague, J. E.: Computing Core-Yeild Limits for Case-Haredened Rolling Bearings,ASLE Trans., 28 (1985) 1, 91-96
3Kunc, R.; Prebil, I.; Torkar, M.:Kovine zlit. tehnol., 34 (2000) 1-2, 31-35
4Lemaitre, J.; A course on damage mechanics,Springer-Verlag, 1996
5Lemaitre, J., Dufailly, J.; Damage measurements, Engineering Fracture Mechanics, 28 (1987) 5/6, 643-661
6Lemaitre, J., Chaboche, J. L.; Mechanics of solid materials, Chambridge University Press, 1990
7Bari, S., Hassan, T.; Anatomy of coupled constitutive models for ratcheting simulation,International Journal of Plasticity,16 (2000), 381-409
8Ohno, N., Wang, J. D.; Kinematic hardening rules with critical state of dynamic recovery, Part I,International Journal of Plasticity, 9 (1993), 375-390
9Ohno, N., Wang, J. D.; Kinematic hardening rules with critical state of dynamic recovery, Part II,International Journal of Plasticity,9 (1993), 391-403
10Armstrong, P. J., Frederick, C. O.; A mathematical representation of the multiaxial Baushinger effect,CEGB Report RD/B/N 731, 1966, Berkeley Nuclear Laboratories
11Petersen, T. O.; Cyclic Plasticity and low cycle fatigue in tool materials, Ph.D. Thesis, Dept. of Solid Mechanics, Technical university of Denmark, DCAMM Report S82, 1998
12Nielsen, H.S.; Intergranular fracture under creep fatigue interaction, Ph.D. Thesis, Dept. of Solid Mechanics, Technical university of Denmark,DCAMM Report S68, 1994
13Calloch, S., Marquis, D.; Triaxial tension-compression tests for multiaxial cyclic plasticity,International Journal of Plasticity, 15 (1999), 521-549
14Abdel-Karim, M., Ohno, N.; Kinematic hardening model suitable for ratchetting with steady-state,International Journal of Plasticity, 16 (2000), 225-240
15Benallal, A., Billardon, R., Lemaitre, J.; Continuum damage mechanics and local approach to fracture: numerical procedures, Computer Methods in Applied Mechanics and Engineering, 92 (1991), 141-155
16Chaboche, J. L.; Continuum damage mechanics: part I - Ge ne ral concepts,Journal of Applied Mechanics, 55 (1988), 59-64
17Chaboche, J. L.; Continuum damage mechanics: part II - damage growth,Journal of Applied Mechanics, 55 (1988), 65-72
18Jiang, Y.; Kurath, P.; Characteristics of the Armstrong-Frederick typeplasticity models,International Journal of Plasticity, 12 (1996) 3, 387-415
19Korelc, J.; Automatic generation of finite-element code by simultaneous optimization of expressions,Theoretical Computer Science,187(1997), 231-248
20Korelc, J.; Automatic derivation of sensivity terms for mixed Galerkin-collocation problems.Zeitschrift für angewandte Mathe- matik und Mechanik, 80 (2000) 2, 533-534
