Numerično modeliranje dodajanja materiala z navarjanjem žice

(1)

UNIVERZA V LJUBLJANI Fakulteta za strojništvo

Numerično modeliranje dodajanja materiala z navarjanjem žice

Magistrsko delo magistrskega študijskega programa II. stopnje Strojništvo

Dejan Kovšca

Ljubljana, avgust 2021

(2)

(3)

(4)

(5)

UNIVERZA V LJUBLJANI Fakulteta za strojništvo

Numerično modeliranje dodajanja materiala z navarjanjem žice

Magistrsko delo magistrskega študijskega programa II. stopnje Strojništvo

Dejan Kovšca

Mentor: izr. prof. dr. Nikolaj Mole Somentor: doc. dr. Bojan Starman

Ljubljana, avgust 2021

(6)

(7)

(8)

(9)

v

Zahvala

Na tem mestu bi se rad zahvalil mentorju izr. prof. dr. Nikolaju Moletu in somentorju doc.

dr. Bojanu Starmanu za strokovno pomoč pri izdelavi te magistrske naloge in vloženo delo pri projektu numeričnega modeliranja dodajanja materiala. Zahvala gre tudi laboratoriju LNMS za izračun numerične simulacije na njihovih računalnikih ter asist. Aljažu Ščetincu in laboratoriju LAVAR za izvedbo eksperimenta oblikovnega navarjanja žice.

Iskreno bi se zahvalil še družini in prijateljem za podporo med študijem.

Želel bi izraziti spomin in posebno zahvalo za življenjski vzor svojemu stricu Jordanu, ki je preminil ravno v času zaključevanja te magistrske naloge.

(10)

vi

(11)

vii

(12)

viii

(13)

ix

Izvleček

UDK 004.942:621.791.75(043.2) Tek. štev.: MAG II/978

Numerično modeliranje dodajanja materiala z navarjanjem žice

Dejan Kovšca

Ključne besede: numerično modeliranje metoda končnih elementov

analiza termo-mehanskega problema dodajanje materiala

dodajne tehnologije obločno navarjanje žice

Tehnološki postopek izgraditve izdelkov po postopku oblikovnega navarjanja žice omogoča izdelavo geometrijsko zahtevnih izdelkov, pri čemer je v primerjavi s konvencionalnimi tehnologijami izkoristek materiala bistveno boljši. Da se določitev procesnih parametrov izdelave izdelkov po tem postopku ne bi opirala samo na metodo poskušanja, se je začel razvoj numeričnih modelov termo-mehanskega odziva izdelka, namenjenih izvedbi računalniške simulacije procesa. Tako lahko ocenimo primernost procesnih parametrov v virtualnem okolju. V okviru te magistrske naloge je obravnavan postopek numeričnega modeliranja oblikovnega dodajanja materiala z obločnim navarjanjem žice na podlagi metode končnih elementov (MKE) v programskem okolju SIMULIA Abaqus. Rezultati računalniške simulacije izdelave geometrijsko enostavnega izdelka so na koncu tudi eksperimentalno ovrednoteni.

(14)

x

(15)

xi

Abstract

UDC 004.942:621.791.75(043.2) No.: MAG II/978

Numerical modelling of material addition by wire welding

Dejan Kovšca

Key words: numerical modelling finite element method thermo-mechanical analysis material addition

additive manufacturing wire arc welding

Wire arc welding-based additive manufacturing (WAAM) is a 3D printing technology for production of near-net-shape parts with complex geometry, where in comparison to conventional production technologies the efficiency of used material is significantly improved. Because the process itself involves thermo-mechanically complex physical phenomena, finite element-based virtual models are commonly employed to optimize the process parameters. This master's thesis presents advanced computational modelling of the WAAM process, based on the SIMULIA Abaqus software environment. By simulating and manufacturing a geometrically simple product, the results of the computer simulation are also experimentally assessed.

(16)

xii

(17)

xiii

Kazalo

Kazalo slik ... xv

Kazalo preglednic ... xvii

Seznam uporabljenih simbolov ... xix

Seznam uporabljenih okrajšav ... xxi

1 Uvod ... 1

1.1 Ozadje problema ... 1

1.2 Cilji ... 1

2 Teoretične osnove in pregled literature ... 3

2.1 Pregled dodajnih tehnologij ... 3

2.2 Metoda končnih elementov ... 4

2.2.1 Toplotni problem ... 4

2.2.1.1 Volumski KE za toplotni problem ... 4

2.2.1.2 Izoparametrični volumski KE za toplotni problem ... 6

2.2.1.3 Gaussova integracijska formulacija ... 8

2.2.1.4 Časovno odvisni toplotni problem ... 9

2.2.2 Mehanski problem ... 10

2.2.2.1 Volumski KE za mehanski problem ... 10

2.2.2.2 Izoparametrični volumski KE za mehanski problem ... 12

2.2.2.3 Selektivna Gaussova integracija mehanskega problema ... 13

2.2.2.4 Iterativno reševanje nelinearnih problemov v plastomehaniki ... 15

2.2.2.5 Konstitutivni model ... 17

2.2.3 Sklopitev termo-mehanskega problema... 25

2.3 Numerično modeliranje dodajanja materiala ... 26

2.3.1 Formulacija računskega problema ... 27

2.3.2 Pristopi modeliranja dodajanja materiala po MKE ... 27

2.3.3 Pogoj za aktivacijo posameznega KE ... 28

2.3.4 Materialne lastnosti neaktiviranih KE ... 29

3 Potek numerične simulacije dodajanja materiala ... 31

3.1 Lastna programska koda iz Pythona/Fortrana ... 33

3.1.1 Določitev aktivacijskih časov ... 33

3.1.2 Iskanje sosednjih KE ... 34

3.1.3 Generacija vhodnih datotek in podprogramov ... 35

(18)

xiv

3.2 Toplotni del izračuna ... 36

3.2.1 Aktivacija končnih elementov s podprogramom USDFLD ... 36

3.2.2 Aktivacija robnih pogojev s podprogramom FILM... 37

3.2.3 Modeliranje toplotnega izvora s podprogramom DFLUX... 38

3.2.4 Poenostavljen grafični prikaz delovanja podprogramov ... 39

3.2.5 Časovno glajenje rezultatov toplotnega izračuna ... 39

3.3 Mehanski del izračuna ... 40

3.3.1 Aktivacija končnih elementov s podprogramom USDLFD ... 41

4 Metodologija raziskave ... 43

4.1 Eksperimentalni del ... 43

4.1.1 Uporabljena oprema za WAAM ... 43

4.1.2 Varilni parametri in opis postopka navarjanja... 44

4.1.3 Meritev temperature pri navarjanju ... 45

4.2 Numerični izračun ... 46

4.2.1 Geometrijski model ... 46

4.2.2 Mreža KE... 47

4.2.3 Določitev aktivacijskih časov z lastno programsko kodo ... 48

4.2.4 Parametri računskih korakov ... 48

4.2.5 Začetni in robni pogoji ... 49

4.2.6 Materialne lastnosti ... 49

5 Rezultati ... 51

5.1 Primerjava rezultatov eksperimenta in numeričnega izračuna ... 51

5.2 Analiza toplotnega stanja ... 52

5.3 Analiza mehanskega stanja ... 54

6 Diskusija ... 57

7 Zaključki ... 59

Literatura ... 61

(19)

xv

Kazalo slik

Slika 2.1: Preslikava izoparametričnega osemvozliščnega heksaedričnega KE iz (a) Kartezijevega v

(b) naravni koordinatni sistem [2] ... 6

Slika 2.2: Primerjava elasto-plastičnosti (a) brez utrjevanja, (b) z linearnim in (c) z nelinearnim utrjevanjem [5] ... 18

Slika 2.3: Primerjava (a) izotropnega, (b) kinematičnega in (c) mešanega utrjevanja [5] ... 19

Slika 2.4: Definicija termičnega raztezka α(T) v ε^α-T diagramu [6] ... 26

Slika 2.5: Primerjava temperaturnega polja med pristopoma dodajanja materiala v Abaqus/Standard: (a) s funkcijo Model Change, (b) z uporabniškim podprogramom USDFLD ... 28

Slika 3.1: Diagram poteka lastne numerične simulacije dodajanja materiala ... 32

Slika 3.2: Primerjava aktivacije sklopa KE v navarku (a) pred naknadno obdelavo aktivacijskih časov in (b) po njej ... 34

Slika 3.3: Geometrija in porazdelitev moči Goldakovega izvora toplote [16] ... 38

Slika 3.4: Grafični prikaz numeričnega modeliranja dodajanja materiala po MKE ... 39

Slika 4.1: Uporabljena oprema pri eksperimentu oblikovnega navarjanja žice ... 44

Slika 4.2: Cilindrična struktura, izdelana s tehnološkim postopkom oblikovnega navarjanja žice . 45 Slika 4.3: Merjenje temperature površine zgornje plasti s termoparom ... 46

Slika 4.4: Geometrijski model obravnavanega cilindričnega izdelka v preseku ... 47

Slika 4.5: Mreža KE obravnavane geometrije v prerezu ... 47

Slika 4.6: Grafični prikaz g-kode, ustvarjene v programskem okolju Ultimaker Cura [14] ... 48

Slika 5.1: Primerjava izračunanih (MKE) in izmerjenih (eksp.) temperatur 20 mm od mesta začetka navarjanja 1., 4., 7. in 10. plasti ... 52

Slika 5.2: Časovni razvoj izračunane temperature v točki 20 mm od mesta začetka navarjanja za posamezno plast ... 53

Slika 5.3: Grafični prikaz časovnega razvoja izračunanega temperaturnega polja pri navarjanju 10. plasti ... 53

Slika 5.4: Grafični prikaz izračunanih radialnih pomikov ur (a) pred odstranitvijo podložne plošče pri sobni temperaturi in (b) po njej ... 54

Slika 5.5: Grafični prikaz izračunanih obodnih napetosti σφφ (a) po končanem navarjanju 10. plasti, (b) po ohladitvi na sobno temperaturo in (c) po odstranitvi podložne plošče ... 55

Slika 5.6: Grafični prikaz časovnega razvoja izračunane Misesove primerjalne napetosti σMises med navarjanjem 10. plasti, po ohlajanju in odstranitvi podložne plošče ... 56

(20)

xvi

(21)

xvii

Kazalo preglednic

Preglednica 2.1: Uteži wi in položaji integracijskih točk ξi Gaussove integracijske formule ... 8

Preglednica 3.1: Prikaz naknadne obdelave aktivacijskih časov na primeru ... 34

Preglednica 3.2: Ploskve in pripadajoča vozlišča heksaedričnega KE v Abaqus/CAE [6] ... 34

Preglednica 4.1: Kemična sestava varilne žice G3Si1 in podložne plošče S235 ... 45

Preglednica 4.2: Meja tečenja pri ekvivalentni plastični deformaciji in temperaturi za material G3Si1 [9] ... 49

Preglednica 4.3: Uporabljene temperaturno odvisne snovne lastnosti za material G3Si1 [9] ... 50

(22)

xviii

(23)

xix

Seznam uporabljenih simbolov

Oznaka Enota Pomen

c J/(kg K) specifična toplota

E N/m² modul elastičnosti

F N sila

f N/m³ volumsko porazdeljena sila

h W/(m² K) prestopnostni koeficient

H N/m² utrjevalni modul

I A električni tok

k W/(m K) koeficient toplotne prevodnosti

P J/s, W moč

q W/m² gostota toplotnega toka

QV W/m³ volumska generacija toplote

T K temperatura

t s čas

U V električna napetost

u m pomik

α m/(m K) razteznostni koeficient

ε / deformacija

𝜀̅^pl / ekvivalentna plastična deformacija

εs / emisivnost

ε^α / termični raztezek

η / izkoristek

ν / Poissonov količnik

ρ kg/m³ gostota

σ N/m² napetost

σMises N/m² Misesova primerjalna napetost

σs W/(m² K⁴) Stefan-Boltzmannova konstanta

σt N/m² meja tečenja

Indeksi

∞ okoliški

dev deviatorični

dil dilatacijski

el elastični

k konvektivni

maks največji

not notranji

pl plastični

pos poskusni

s sevalni

zun zunanji

(24)

xx

(25)

xxi

Seznam uporabljenih okrajšav

Okrajšava Pomen

3D tridimenzionalno

CAD računalniško podprto načrtovanje (angl. computer-aided design) CAM računalniško podprta proizvodnja (angl. computer-aided

manufacturing)

CNC računalniško numerično krmiljenje (angl. computerized numerical control)

KE končni element

KKS Kartezijev koordinatni sistem KTO konsistentni tangentni operator

MIG varjenje kovin z uporabo inertnih plinov (angl. metal inert gas)

MKE metoda končnih elementov

NKS naravni koordinatni sistem

WAAM oblikovno navarjanje žice (angl. wire arc additive manufacturing)

(26)

xxii

(27)

1

1 Uvod

1.1 Ozadje problema

Dodajanje materiala oz. t. i. 3D-tisk uvrščamo med dodajne tehnologije. Gre za skupino tehnoloških postopkov, za katere je značilno oblikovno dodajanje materiala po plasteh. Ti postopki se uveljavljajo predvsem na področju izdelave kompleksnih geometrij, ki jih s konvencionalnimi postopki ne moremo izdelati. Dodatno njihova uporaba omogoča bistveno boljši izkoristek materiala.

Razvoj numeričnih modelov termo-mehanskega odziva izdelka pri dodajanju materiala na podlagi metode končnih elementov (MKE) omogoča določitev in oceno primernosti procesnih parametrov v virtualnem okolju. Na ta način se lahko v praksi izognemo stroškovno neupravičeni metodi poskušanja (angl. trial and error approach). Poleg tega pa nam računalniška simulacija omogoča napoved deformacijskega in napetostnega stanja v izdelku po koncu izdelave.

1.2 Cilji

V okviru te naloge bodo najprej opisane teoretične osnove preračuna termo-mehanskega problema po MKE in splošni pristopi k numeričnemu modeliranju dodajanja materiala. V nadaljevanju bo predstavljen lasten numerični model dodajanja materiala z obločnim navarjanjem žice v programskem okolju SIMULIA Abaqus. Na podlagi razvitega modela bosta izvedeni računalniška simulacija izdelave geometrijsko enostavnega izdelka in analiza termo-mehanskega stanja v njem. Rezultati bodo na koncu tudi eksperimentalno ovrednoteni.

(28)

Uvod

2

(29)

3

2 Teoretične osnove in pregled literature

2.1 Pregled dodajnih tehnologij

V splošnem računalniško vodeni premični izvor toplote po plasteh pretaljuje dodajani material. Najpogosteje se za proizvodnjo uporabljajo polimeri ali kovine. Glede na način dodajanja kovinskih materialov v inertni atmosferi ločimo tri glavne postopke [1]:

‐ Selektivno pretaljevanje prahu: Dozirnik enakomerno razdeli kovinski prah po celotni površini delovne mize (angl. powder bed). Laserski ali elektronski žarek pa nato selektivno pretali kovinski prah zgolj na predvidenih mestih. Ko je taljenje trenutne plasti zaključeno, se miza spusti za debelino nove plasti. Celoten postopek se ponavlja, dokler izdelek ni končan. Prednost opisanega postopka je visoka dimenzijska natančnost, po drugi strani pa je ta postopek zato relativno počasen in drag. Velikost izdelka je omejena z velikostjo delovne mize.

‐ Direktna depozicija žice: Gre za oblikovno navarjanje varilne žice, pri katerem je vir dodane toplote lahko energija varilnega obloka, laserski ali elektronski žarek. Od omenjenih postopkov ta omogoča največjo hitrost dodajanja materiala in je posledično manj natančen. Zaradi naštetega se uporablja za večje izdelke in zahteva večjo stopnjo nadaljnjih obdelav. Hkrati pa zahteva najnižjo investicijo v opremo in omogoča uporabo cenovno najugodnejšega materiala.

‐ Direktna depozicija prahu: Laserski ali elektronski žarek na površini ustvarja talino prahu, ki se skozi šobe dodaja neposredno na mesto navarjanja. Cena, dimenzijska natančnosti in hitrost dodajanja materiala ga uvrščajo med oba zgoraj opisana postopka.

Direktna depozicija prahu med drugim omogoča tudi navarjanje na že obstoječe izdelke in prehod med različnimi dodajnimi materiali.

Pri vseh omenjenih postopkih je zaradi zahtevnosti dodajanja kovinskih materialov potrebna naknadna obdelava. Ta zajema konvencionalne procese odrezovanja, kot so frezanje, brušenje ali poliranje z abrazivnim tokom. Z njimi vplivamo na dimenzijsko natančnost ter stanje (hrapavost) in videz zunanjih površin. Dodatno lahko dosežemo tudi zmanjšanje površinskih zaostalih napetosti.

(30)

Teoretične osnove in pregled literature

4

2.2 Metoda končnih elementov

Numerične simulacije termo-mehanskih problemov v praksi pogosto temeljijo na metodi končnih elementov (MKE). Ta je zasnovana na šibki obliki integralske formulacije problema. Prednosti slednje sta izpolnjevanje vodilne diferencialne enačbe v povprečju po celotnem območju in znižan red diferencialnega operatorja vodilne enačbe problema.

MKE aproksimira vrednosti polja primarne spremenljivke z neznanimi vrednostmi v vozliščih končnih elementov (KE). Pri reševanju problema moramo poiskati rešitev sistema toliko enačb, kolikor ima problem prostostnih stopenj. Število prostostnih stopenj problema je enako produktu vozlišč, ki smo jih pri diskretizaciji geometrijskega območja definirali, ter številu neznank v posameznem vozlišču. Izoparametrični KE nam pri tem omogočajo obravnavo geometrijsko zahtevnejših problemov.

2.2.1 Toplotni problem

V tem poglavju bo predstavljena uporaba MKE v analizi toplotnega problema. Izpeljali bomo enačbe osemvozliščnega izoparametričnega volumskega KE za stacionarni in časovno odvisni problem prevoda toplote.

2.2.1.1 Volumski KE za toplotni problem

Izhodišče za izpeljavo enačbe splošnega volumskega KE za toplotni problem [2] je vodilna enačba stacionarnega prevoda toplote (2.1), kjer T [K] predstavlja temperaturno polje, k [W/(m K)] toplotno prevodnost in QV [W/m³] volumsko generacijo toplote. Zaradi preglednosti je uporabljen tenzorski zapis enačb, kjer je parcialni odvod po j-ti komponenti označen z indeksom », j«:

( )

,_j , V 0 , , ,

k T j+Q = j=x y z . _(2.1)

Vodilno parcialno diferencialno enačbo pomnožimo s poljubno skalarno funkcijo v in produkt integriramo po območju Ω posameznega KE. Poljubna funkcija v nam omogoča, da dobimo potrebno število medsebojno neodvisnih enačb. Za enolično rešitev namreč potrebujemo enako število enačb in neznank. Po preureditvi dobimo:

(

,_j

)

, d V d 0 , , ,

k T jv Q v j x y z

 

 +  = =

 

^. (2.2)

Ob upoštevanju matematične zveze za parcialno odvajanje sestavljene funkcije (2.3):

(

,_j

) ( )

, ,_j , ,_j ,_j , , ,

j j

k T v = k T v k T v+ j=x y z , _(2.3)

(31)

5 lahko enačbo (2.2) zapišemo kot:

( )

,_j ,_j d ,_j , d V d , , ,

k T v k T v j Q v j x y z

  

   =  +  =

 

  

^. ^(2.4)

Z uporabo Gaussovega divergenčnega teorema (2.5) znižamo red diferencialnega operatorja in integral po volumnu Ω prevedemo na integral po površini območja Γ, ki ga omejuje. Pri čemer nj označuje komponente enotskega normalnega vektorja na površino Γ:

(

,_j

)

, d ,_j _j d _j _j d , , ,

j v

I k T v k T n q n v j x y z

  

   

=



 =



   = −



   = ^. ^(2.5)

Upoštevajoč Fourierjev zakon za gostoto toplotnega toka qj = -k T,j za j = x,y,z, izpeljemo šibko obliko integralske formulacije (2.6), ki je osnova za MKE:

,_j ,_j d _j _j d V d , , ,

k T v q n v Q v j x y z

  

   = −    +  =

   

  

^. ^(2.6)

V skladu z Galerkinovo metodo izberemo poljubne funkcije v enake interpolacijskim funkcijam ψ, in zapišemo:

,_j _{i j}, d _j _j _i d V _id , , , , 1,..., v

k T q n  Q  j x y z i N

  

   = −    +  = =

   

  

^. ^(2.7)

Temperaturno polje T(x,y,z) znotraj območja posameznega KE aproksimiramo z interpolacijo prek vozliščnih vrednosti temperature Ti (2.8). V ta namen uporabimo polinomske interpolacijske funkcije ψ. Pri čemer indeks i označuje posamezno funkcijo, ki pripada določenemu vozlišču KE, tako da imamo skupaj Nv različnih funkcij. Z {N}

označimo vektor interpolacijskih funkcij, s {T} pa vektor vozliščnih temperatur:

( ) ( )

^v

( )    

^T

1

, , ˆ , , , ,

N

i i

i

T x y z T x y z T  x y z N T

=

 =



 = ^. ^(2.8)

Za nadaljnjo lažjo izpeljavo uvedemo vektor odvodov {L} = {^∂/∂x,^∂/∂y,^∂/∂z}^T in zgornjo integralsko enačbo zapišemo v vektorski obliki (2.9), kjer {q} predstavlja vektor gostote toplotnega toka normalno na površino Γ in {Q} vektor volumske generacije toplote v območju Ω:

  

( )

^T

  

^d

     

^T ^d

   

^T ^d

k L N L N T N q N Q

  

 

 = −  + 

 







 

^. ^(2.9)

(32)

6

2.2.1.2 Izoparametrični volumski KE za toplotni problem

Zato da topologija mreže KE v splošnem omogoča ustreznejši popis obravnavane geometrije in da je lažje implementirati Gaussovo numerično integracijo, vpeljemo izoparametrične KE. Pri tem heksaedrične KE poljubnih oblik v Kartezijevem koordinatnem sistemu (KKS) preslikamo v naravni koordinatni sistem (NKS), kjer vsi zavzemajo obliko kocke s središčem v točki (0,0,0) in stranico dolžine dveh enot (glej sliko 2.1).

(a) (b)

Slika 2.1: Preslikava izoparametričnega osemvozliščnega heksaedričnega KE iz (a) Kartezijevega v (b) naravni koordinatni sistem [2]

Za interpolacijo geometrije in primarnih spremenljivk (temperatura, pomiki v posamezni koordinatni osi) po območju KE uporabimo polinomske interpolacijske funkcije ψ, ki so za volumski osemvozliščni heksaedrični KE v NKS (ξ,η,μ) definirane z enačbo (2.10).

Posamezno vozlišče KE označuje indeks j:

( ^{, ,} )

(

¹

)(

¹

)(

¹

)

^{/ ,}⁸ ¹^,...,⁸

j j j j j

    = +   +   +   = . (2.10)

V primeru izoparametričnega KE geometrijo v KKS (x,y,z) interpoliramo z vozliščnimi vrednostmi Xⁱ v odvisnosti od interpolacijskih funkcij  v NKS (ξ,η,μ) na naslednji način:

( ) ( )

1 Nv

i

i i j j

j

x x    X     i x y z

=

= ^{, ,} 



 ^{, , ,} = ^{, ,} ^. ^(2.11)

Ker izvajamo preslikavo iz KKS v NKS, nas bodo v nadaljevanju zanimali parcialni odvodi Kartezijevih koordinat (x,y,z) po naravnih koordinatah (ξ,η,μ), ki jih zapišemo z Jacobijevo matriko [J]:

 

x x x

y y y

z z z

  

  

    

 

  

 

=    

   

  

 

 

J . (2.12)

(33)

7 Povezavo med volumnom KE v KKS in NKS predstavlja determinanta Jacobijeve matrike:

( )  

d d

 

 = 

 

^det ^J ^. ^(2.13)

Zvezo med delom površine KE v KKS in NKS pa predstavlja determinanta dela Jacobijeve matrike (2.14). Ta je za primer stranice z normalo v smeri μ izračunana z (2.15):

( )  

d d

 

 = 

 

^det ^j ^, ^(2.14)

e e e

x y z

y y y

  

   

  

 

 

   

  =

      

   

    

 

j . _(2.15)

Relacijo med parcialnimi odvodi interpolacijskih funkcij v KKS {N} po kartezijevih koordinatah in parcialnimi odvodi interpolacijskih funkcij v NKS {Ñ} po naravnih koordinatah izrazimo prek inverza Jacobijeve matrike [J]^-1. Parcialne odvode interpolacijskih funkcij krajše zapišemo z matriko [B] (2.16):

 

^B ⁼  ^{L N} ⁼

 

^J ^-1 ^{L N}

 

^. ^(2.16)

Končno lahko vse zgoraj izpeljane zveze vstavimo v enačbo (2.9) in zapišemo:

   

^T ^det

( )  

^d  

 

^T ^det

( )  

^d

 

^T ^de^t

( )  

^d

k T N q N Q



 

 

 = −  + 

 





^B ^B ^J 



^j



^J ^, ^(2.17)

oziroma krajše:

 

M e

 

T e =

   

qq e + qQ e . (2.18) Kjer [M]e predstavlja prevodnostno matriko posameznega KE, {T}e vektor pripadajočih vozliščnih temperatur, {qq}e vektor ekvivalentnih vozliščnih izvorov/ponorov toplote zaradi prevoda toplote preko roba KE in {qQ}e vektor ekvivalentih vozliščnih izvorov/ponorov zaradi volumske generacije toplote po območju KE.

Za posamezni KE dobimo toliko enačb, kolikor ima KE prostostnih stopenj. V primeru volumskega KE za toplotni problem je v vsakem vozlišču zgolj ena neznana primarna veličina – neznana vozliščna temperatura Ti. Skupno ima tako posamezni KE Nv

prostostnih stopenj. Če sistem linearnih enačb razširimo na vse KE, ki popisujejo

(34)

8

obravnavano območje problema, moramo komponente togostne matrike posameznih KE sešteti. Tako dobimo globalno prevodnostno matriko [M]:

 

M

 

T =

   

qq + qQ . (2.19) Zgornji sistem linearnih enačb lahko rešujemo direktno ali iterativno [6]. Direktno reševanje z Gaussovo eliminacijo je sicer računsko intenzivnejše, vendar praviloma daje numerično natančnejšo rešitev. Hkrati je ta, tudi iz vidika pričakovanega števila neznank v obravnavanem primeru, ustreznejša izbira.

2.2.1.3 Gaussova integracijska formulacija

Pri MKE integrale v enačbi (2.17) praviloma računamo z Gaussovo integracijsko formulo (angl. quadrature rule). Ta numerična metoda integriranja je bistveno hitrejša od analitičnega integriranja in preostalih numeričnih metod, saj temelji na izračunu vsote uteženih vrednosti funkcije v diskretnih točkah (2.20). Vrednosti uteži wi in naravne koordinate ξi so za izbrano število integracijskih točk m prikazane v preglednici 2.1. Ker formulacija temelji na polinomski aproksimaciji, pri MKE pa uporabljamo polinomske interpolacijske funkcije, lahko z m integracijskimi točkami v eni dimenziji izračunamo eksaktno integral polinoma n = 2m - 1 stopnje:

( )

1

1 1

d

m

i i

i

f x x w f 

+

− =









^. ^(2.20)

Preglednica 2.1: Uteži wi in položaji integracijskih točk ξi Gaussove integracijske formule

m i wi ξi

1 1 2 0

2 1 1 ¹ ₃

2 1 - ¹ ₃

Če v trirazsežnem prostoru uporabimo dve integracijski točki v posamezni dimenziji, dobimo za heksaedrični osemvozliščni KE osem integracijskih točk, v primeru reducirane integracije pa zgolj eno. Numerično integracijo v prostoru izvedemo kot vsoto uteženih vsot v posamezni koordinatni smeri.

Prevodnostno matriko [M]e in vektorja ekvivalentnih vozliščnih izvorov {qq}e in {qQ}e za posamezni KE z Gaussovo integracijo izračunamo na podlagi vsote j = 1, …, 8 integracijskih točk (m = 2), z utežjo w = 1, po naslednjih enačbah:

  ⁸    ^{e T} ^e

( )

 

1

e

e e j j

j

k

=



M B B det J , _(2.21)

(35)

9

 

⁴

 

^T^{ }

( )

 

1

e e

q e j

j

q N q

=



^det ^j ^, ^(2.22)

 

⁸

 

^T^{ }

( )

^{ }

1

e e

Q e j

j

q N Q

=



^det ^J ^. ^(2.23)

2.2.1.4 Časovno odvisni toplotni problem

V primeru nestacionarnega toplotnega problema moramo upoštevati časovno odvisno vodilno enačbo prevoda toplote (2.24). Ta ima, za razliko od enačbe stacionarnega prevoda toplote (2.1), dodaten, časovno odvisni člen (parcialni odvod temperature po času): ρcT,t. Vsebuje pa še dve snovni lastnosti: gostoto ρ [kg/m³] in specifično toploto c [J/(kg K)]:

(

^{k T}^,^j

)

_,_j⁺^Q^V ⁼^^cT^,^t ^, ^j⁼^{x y z}^{, ,} ^. (2.24)

Izpeljava enačbe KE za nestacionarni prevod toplote (2.26) je ekvivalentna primeru stacionarnega problema (glej poglavje 2.2.1.1 Volumski KE za toplotni problem). Pri tem pa parcialni odvod temperature aproksimiramo z diferenčno shemo (2.25):

( , , , ) ( , , , _k 1) ( , , , _k)

T x y z t T x y z t T x y z t

t t

 + −

   , (2.25)

 

, , V

v

1

d d d d ,

, , , 1,..., , 0,1

.

k k k

j I j j j I I

k k

k T q n Q

t j x y z I

T

N

cT

    



+ + +

   

   = −    +  − + 

    

=

−

= 

   

_(2.26)

Z izbiro koeficienta β ∈ [0,1] določamo, v katerem časovnem trenutku bomo diferencialno enačbo upoštevali. V enačbi KE uporabimo aproksimacijo veličin v časovnem trenutku tk+β, medtem ko je neznana veličina, ki jo računamo, vrednost temperature v časovnem trenutku tk+1. Diskretne vrednosti spremenljivk T^k, q^k, Q^kV = f(x,y,z,tk) v časovnem trenutku tk+β določimo z linearno interpolacijo vrednosti pri času tk in tk+1 (2.27):

(

, , , _k

)

⁽ , , , _k⁾⁽1 ⁾ ⁽ , , , _k 1⁾

f x y z t ₊_ = f x y z t − + f x y z t ₊  . _(2.27)

Glede na izbiro vrednosti koeficienta β ločimo tri metode reševanja časovno odvisnih problemov:

‐ eksplicitna metoda: β = 0,

‐ implicitna metoda: β = 1,

‐ Crank-Nicolsonova metoda: β = 0,5.

(36)

10

2.2.2 Mehanski problem

V tem poglavju bo predstavljena uporaba MKE v kvazistatičnih nelinearnih mehanskih analizah. Izpeljali bomo enačbe osemvozliščnih izoparametričnih volumskih KE in jih skupaj z ustreznimi konstitutivnimi enačbami plastičnosti implementirali v algoritem za izračun elastoplastičnega problema.

2.2.2.1 Volumski KE za mehanski problem

Izpeljava sistema enačb splošnega volumskega KE za mehanski problem [2] sledi iz parcialnih diferencialnih enačb ravnotežnega stanja v posamezni smeri (2.28), kjer σij [MPa] predstavljajo komponente napetostnega tenzorja, fi [N/m³] pa volumsko porazdeljeno obremenitev (npr. lastna teža). Izračunamo jo kot produkt gostote ρ in pospeška ai, po enačbi: fi = ρ ai. Zaradi preglednosti je uporabljen tenzorski zapis enačb:

ij j fi 0 i x y z

 _, + = , = , , . _(2.28)

Zgoraj zapisane parcialne diferencialne enačbe pomnožimo s poljubnimi funkcijami vi in integriramo po območju Ω posameznega KE. Po preureditvi dobimo:

d d 0

ij j iv f vi i i x y z



 

 +  = =



^,



^, ^{, ,} ^. ^(2.29)

Ob upoštevanju matematične zveze za parcialno odvajanje sestavljene funkcije (2.30):

(

_{ij i}

)

, _{ij j i}, _{ij i j}, , , ,

v j v v i x y z

 = + = , _(2.30)

lahko enačbo (2.29) zapišemo kot:

d d d

, ( ), , , ,

ij i jv ij iv j f vi i i x y z

 

  

 =  +  =

  

^. ^(2.31)

Z uporabo Gaussovega divergenčnega teorema znižamo red diferencialnega operatorja ter enačbo (2.31) prevedemo na vsoto integrala po volumnu Ω in integrala po površini območja Γ, ki ga omejuje:

( )

d d d

ij i jv ij jn vi f vi i i x y z

 

  

 =  +  =



^,

 

^, ^{, ,} ^. ^(2.32)

Upoštevajoč definicijo napetostnega vektorja pi = σij nj, zapišemo:

d d d

ij i jv p vi i f vi i i x y z



  

 =  +  =



^,

 

^, ^{, ,} ^. ^(2.33)

(37)

11 Zgornji tenzorski zapis (2.33) predstavlja tri integralske enačbe, ki jih lahko medsebojno seštejemo, saj zajemajo isto območje Ω in površino Γ, ki to območje ograjuje. Tako nastali sistem enačb vektorsko zapišemo:

  

^ ^v ^d

   

^p ^T ^v ^d

   

^f ^T ^v ^d

  

  =  + 

  

^. ^(2.34)

Zaradi simetričnosti uporabimo Voigtov zapis napetostnega σij in deformacijskega εij

tenzorja s šestimi komponentami:

 ^ =



     ^xx^, ^yy^, ^zz^, ^xy^, ^xz^, ^yz



^T^, ^{ } =



     ^xx^, ^yy^, ^zz^, ^xy^, ^xz^, ^yz



^T . (2.35)

Zvezo med napetostnim in deformacijskim tenzorjem v splošnem zapišemo z enačbo (2.36). Matrika [D] predstavlja t. i. materialno matriko oziroma konsistentni tangentni operator (glej poglavje 2.2.2.5 Konstitutivni model). Relacija med deformacijami {ε} in vektorjem pomikov {u} = {ux, uy, uz}^T pa je podana z enačbo (2.37). Pri tem smo definirali matriko odvodov [L] in kompenzacijsko matriko [W], saj za strižne deformacije γij velja:

εij = ½γij, za vsak i ≠ j:

  =

  

D W     =  ^T

  

D W , (2.36)

 

^W

 

^ ⁼

 

^L

 

^u ^, ^(2.37)

 

0 0

0 0 0

x y

z

y x

z x

z y



 

 

 

=  

 

 

L ,

 

1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 2 0 0 0 0 0 0 2

 

 

=  

 

 

W . (2.38)

V primeru homogenega, izotropnega in linearno elastičnega materiala zvezo med napetostmi in deformacijami podaja posplošen Hookov reološki zakon, kjer je [D] podan z enačbo (2.89).

Polje pomikov {u} aproksimiramo z interpolacijo vozliščnih vrednosti pomikov Uⁱ znotraj območja posameznega KE. V ta namen uporabimo interpolacijske funkcije ψ. Pri tem indeks j označuje posamezno od skupno Nv vozlišč KE. Z [N] pa označimo matriko interpolacijskih funkcij:

( ) ( )  

^T

 

1 Nv

i i

i j j

j

u x y z U  x y z  U i x y z

=





 = =

, , , , , , , , (2.39)

(38)

12

 

   

  ^{  }

T T

T

0 0

x x

y y

z z

u U

u u U U

u U



   

     

    

=  =  =

 

   

 

    

N . _(2.40)

Poljubne funkcije v izberemo v skladu z Galerkinovo metodo, kot interpolacijo poljubnih vozliščnih vrednosti ϒ (2.41):

 

   

  ^{  }

T T

T

0 0

x x

y y

z z

v

v v

v



    

     

    

=  =   = 

 

     

    

N . _(2.41)

Zdaj lahko vse zgoraj zapisane zveze vstavimo v enačbo (2.34), izpostavimo poljubne vozliščne vrednosti ϒ in zapišemo:

    

(

^U

)

^T

      ( )

^d

   

^p ^T

 

^d

 

^f ^T

 

^d

 

  

   

   =   +  

   

   





^{L N} ^{D W} ^{L N}  



^N



^N  ^. ^(2.42)

Tako dobljena enačba mora veljati za poljubne vozliščne vrednosti ϒ, kar je v splošnem mogoče le, če velja enakost (2.43). Pri tem smo iz integrala izpostavili konstante vozliščnih vrednosti pomikov {U}:

  

( )

^T

      ( )

^d

^{ }

^U

^{ }

^p ^T

 

^d

 

^f ^T

 

^d

  

 

 =  + 

 





^{L N} ^{D W} ^{L N} 



^N



^N ^. ^(2.43)

2.2.2.2 Izoparametrični volumski KE za mehanski problem

Relacijo med parcialnimi odvodi interpolacijskih funkcij v KKS [N] po kartezijevih koordinatah in parcialnimi odvodi interpolacijskih funkcij v NKS [Ñ] po naravnih koordinatah izrazimo prek inverza Jacobijeve matrike [J]^-1. Parcialne odvode interpolacijskih funkcij krajše zapišemo z matriko [B]:

        

^B ⁼ ^{L N} ⁼ ^J ^-1 ^{L N}^{ }_{ }^. (2.44) Vse že izpeljane zveze v poglavju 2.2.1.2 Izoparametrični volumski KE za toplotni problem, vključno z (2.44), vstavimo v levi del enačbe KE (2.43) in zapišemo:

(39)

13

    

^T

( )   ( )  

^d

 

^U

 

^p ^T

 

^d

 

^f ^T

 

^d

 



 

 =  + 

 





^B ^{D W} ^B ^det ^J 



^N



^N ^, ^(2.45)

oziroma krajše:

 K e U e =Fzune . _(2.46)

V enačbi (2.46) [K]e predstavlja togostno matriko posameznega KE, {U}e vektor pripadajočih vozliščnih pomikov in {Fzun}e vektor ekvivalentih vozliščnih oziroma točkovnih zunanjih sil.

V primeru volumskega KE so v vsakem vozlišču tri neznane primarne veličine – tri komponente vektorja pomika. Skupno ima tako posamezni KE 3Nv prostostnih stopenj. Če sistem enačb razširimo za poljubno število KE, moramo komponente togostne matrike posameznega KE prišteti na ustrezno mesto v globalni togostni matriki [K]:

 

K

 

U =

 

Fzun . (2.47)

Iz zgornje enačbe sledi ravnotežna enačbo sil (2.48), ki enači zunanje in notranje sile v oziroma na posameznem KE. Notranje sile {Fnot}e lahko za vsak posamezni KE iz napetosti izračunamo po enačbi (2.49). Vektor zunanjih sil {Fzun} določimo na podlagi robnih pogojev, bodisi kot znano (prosta vozlišča) ali neznano vrednost (podprta vozlišča):

   

Fnot = Fzun , (2.48)

 

^not

 

^T

  ( )  

^d

e

e e e e e

F 



=



^B ^det ^J  ^. (2.49)

2.2.2.3 Selektivna Gaussova integracija mehanskega problema

Togostno matriko posameznega KE [K]e, ki sledi iz enačbe (2.45), izračunamo po enačbi (2.50). Pri integraciji slednje po območju posameznega KE ugotovimo, da se po prostoru spreminja le vrednost matrike odvodov interpolacijskih funkcij v NKS [B]e. Matrika [D]e, ki popisuje zvezo med napetostmi in deformacijami (t. i. materialni modul), je določena za vsako integracijsko točko posebej, vendar je odvisna zgolj od materialnih lastnosti in je zato ne integriramo po volumnu. Jacobijeva matrika [J]e je za vsak posamezni KE enolično določena, njena vrednost pa se po volumnu KE prav tako ne spreminja. Kompenzacijska matrika [W] pa je konstanta: