MNM Numerične metode modeliranja

Učinek = vzrok, je neodvisna veličina. Odziv = posledica, je od učinka odvisna veličina. Primer: na vzmet obesimo maso $m$ (učinek), raztezek $x = c_1 m$ je odziv. Vlogi se lahko tudi zamenjata.

• Neposredni učinek je učinek, brez katerega ni odziva (npr. sila na vzmet).
• Posredni učinek so ostali parametri sistema, ki tudi vplivajo na odziv (npr. vzmetna togost $k$).

$$\text{odziv} = f(\text{učinek},\ \text{konstitutivni parametri})$$

Predpogoj za uspešno modeliranje je fizikalno razumevanje problema. Odgovoriti moramo na vprašanja:

• Zakaj? npr. Zakaj pridržujemo pločevino pri globokem vleku?
• Kako? npr. Kako določimo vrstni red izrezov lamel v šestrednem orodju?
• Kdaj? npr. Kdaj pride do pretrga pločevine pri rezanju?
• Kje? npr. Kje se ulitek počasneje ohlaja?

Do modela pridemo z opazovanjem pojavov in z ugotavljanjem soodvisnosti med veličinami, ki pojav opredeljujejo (učinek – odziv).

Pojav opazujemo in analiziramo na primerih, ki so za opazovani pojav dovolj signifikantni. Tak primer opredeljuje fizikalni (eksperimentalni) model. Fizikalni model definira območje in čas opazovanja ter objekte, pomembne za objektivno identifikacijo opazovanega pojava.

Razvoj veličin v fizikalnem modelu opredeljujejo naravne zakonitosti: splošni aksiomi (akcija = reakcija), fizikalni zakoni (Newtonovi zakoni) in konstitutivni zakoni (trdnina, kapljevina).

Vzpostavitev odnosov, ki obvladujejo opazovani pojav, opredeljuje matematični model. Sestavljajo ga algebrajske, diferencialne ali integralske enačbe. Vsak matematični model je postavljen v okvir prostora in časa.

Numerični model je računsko obvladovanje matematičnega modela. Enačbe lahko rešujemo na 2 načina:

1. Eksaktno analitično reševanje s funkcijskimi rešitvami v zaključeni obliki. Napake izhajajo iz nenatančnih vhodnih podatkov in iz zapisa realnih števil v računalniški spomin.
2. Aproksimativno numerično reševanje z rešitvami v diskretni obliki (v posameznih točkah območja). Napake izhajajo iz nenatančnih vhodnih podatkov, zapisa realnih števil v računalniški spomin in iz izbrane numerične metode. Med diskretnimi točkami rešitev ekstrapoliramo (linearno, kvadratično).

Eksaktno reševanje da funkcijsko rešitev v zaključeni obliki. Ta zadošča vsem točkam območja in robnim pogojem.

Aproksimativno reševanje da rešitev v diskretnih točkah, med katerimi rešitev ekstrapoliramo. Dodatna napaka izhaja iz izbrane numerične metode.

Napake izhajajo iz:

• nenatančnih vhodnih podatkov,
• zapisa realnih števil v računalniški spomin,
• izbrane numerične metode reševanja.

Z osvojenim numeričnim modelom je omogočeno računalniško simuliranje. Na simulaciji sta zasnovana načrtovanje in odločanje (o procesu, produktu ...).

Simuliranje je iskanje odziva pri enem ali več naborih spreminjajočih se vhodnih podatkov. Iz odzivov za različne učinke s sintezo pridemo do odločitve. Numerični model omogoča tudi avtomatsko krmiljenje tehnoloških procesov.

Glede na gibalno stanje materialnih delcev lahko dogodke obravnavamo v geometrijskem ali materialnem prostoru.

• Geometrijski prostor: opazovanje je vezano na časovno fiksno območje geometrijskih točk v prostoru. Ni pomembno, ali so te točke v različnih trenutkih zasedene z različnimi materialnimi delci (npr. vlečenje žice).
• Materialni prostor: opazovanje je vezano na fiksno območje snovnih točk, ne glede na njihovo gibalno stanje (npr. preoblikovanje v testastem stanju).

Če snovne točke s časom ne spreminjajo položaja, sta geometrijski in materialni prostor identična.

Čas je fizikalna danost:

• čas ima razsežnost (je enorazsežen),
• čas je progresivna spremenljivka, njegova vrednost lahko le narašča.

Trorazsežni materialni prostor in čas tvorita štirirazsežni hiperprostor.

Koordinatni sistem omogoča matematični popis materialnega prostora in dogodkov v njem. Izbira K.S. ne vpliva na fizikalno vsebino problema. Ima pa lahko ključen vpliv na matematično formulacijo in numerično reševanje problema.

Modelirano območje $\Omega^{M}$ je območje našega zanimanja v materialnem prostoru $\Omega^{\infty}$. Običajno je omejeno: $$\Omega^{M} \subset \Omega^{\infty}$$ Ostalo območje je okolica modeliranega območja: $$\Omega^{\infty - M} = \Omega^{\infty} - \Omega^{M}$$ Zaprto podobmočje $\Omega^{M}$ in odprto podobmočje $\Omega^{\infty - M}$ izkazujeta lastnosti: $$\Omega^{M} \cup \Omega^{\infty - M} = \Omega^{\infty} \quad in \quad \Omega^{M} \cap \Omega^{\infty - M} = \emptyset$$ Mejo med njima tvori ograja $\Gamma^{M}$ (robne točke $\Omega^{M}$). Določa jo Dedekindov presek odprtih podobmočij: $$(\Omega^{M} - \Gamma^{M})/\Omega^{\infty - M} = \Gamma^{M}$$ Modelirano območje je lahko sestavljeno iz več zaprtih podobmočij $\Omega^{M} = \bigcup_{k=1}^{N} \Omega_k$. Mejo med sosednjima podobmočjema določa presek $\Omega_i \cap \Omega_j = \Gamma_{ij}$.

Glede na fizikalne in geometrijske posebnosti problema je matematični model velikokrat mogoče zasnovati v geometrijskem prostoru, katerega razsežnost je manjša od razsežnosti materialnega prostora: $$\Re^3 \rightarrow \Re^p \ ; \quad p \in \{1, 2\}$$ Tako dobimo 3D, 2D ali 1D model. Sama geometrija za redukcijo ni dovolj. Redukcijo morajo dopuščati tudi obremenitve in robni pogoji. Redukcija bistveno skrajša čas računanja.

Diskretni sistemi so sistemi s končnim številom prostostnih stopenj.

V modeliranem območju so podobmočja:

• ki ne mejijo druga z drugo,
• katerih medsebojna razdalja je praviloma veliko večja od najdaljše dolžine kateregakoli izmed teh podobmočij,
• katerih porazdelitev snovi nima odločujočega vpliva na odziv preostalih podobmočij; vpliv je odvisen le od količine snovi v podobmočju,
• katerih lastnosti snovi se bistveno razlikujejo od snovi v sosednjih podobmočjih.

Primeri: sistem masnih točk, polja točkovno porazdeljenih električnih nabojev, kristalne kali v talini.

Kontinualni sistemi so sistemi z neskončnim številom prostostnih stopenj.

V modeliranem območju je:

• snov porazdeljena zvezno po posameznih podobmočjih, z morebitno nezveznostjo le na prehodu med sosednjimi podobmočji,
• porazdelitev snovi po podobmočju odločilna za odziv v preostalih podobmočjih.

Kontinualni sistem je limita večprostostnega diskretnega sistema, katerega število prostosti preseže vse meje, razdalje med delci pa so infinitezimalne. Spremenljivke so zvezne funkcije prostorskih koordinat in časa. Reševanje se prevede na reševanje ene ali več diferencialnih enačb.

Primeri:

• plastovite kompozitne plošče
• območje zrak-morje-kopno
• dvofazno območje led-voda (taljenje, strjevanje)
• dvofazno območje trdnina-kapljevina (vodne turbine)
• tehnološko postrojenje za litje

	Diskretni sistemi	Kontinualni sistemi
Število prostostnih stopenj	končno	neskončno
Porazdelitev snovi	ločena podobmočja; na ostale vpliva le količina snovi	zvezna; porazdelitev je odločilna za odziv ostalih podobmočij
Reševanje	algebrajske enačbe ali NDE po času	diferencialne enačbe po prostoru (in času)

Čas je v absolutnem pogledu progresivna spremenljivka, $dt > 0$. Kljub temu lahko nekatere probleme obravnavamo kot časovno nespremenljive.

• Stacionarni modeli (časovno neodvisni): $\frac{d}{dt}(\text{učinek},\ \text{konstitutivni parametri},\ \text{odziv}) = 0$
• Nestacionarni modeli (časovno odvisni): $\frac{d}{dt}(\text{učinek},\ \text{konstitutivni parametri},\ \text{odziv}) \neq 0$

Časovna odvisnost je relativna. Transportni problemi, ki so v krajevnem smislu stacionarni (vlečenje žice, kontinuirno litje), so v materialnem prostoru nestacionarni. Obravnava v geometrijskem prostoru omogoča izločitev časovne dimenzije.

Matematični model naj bo kolikor je mogoče enostaven. Njegova zahtevnost naj bo ravno tolikšna, da zaobjame vse ključne dejavnike za verodostojnost sistemskega odziva.

Odvečna zahtevnost ne da nujno bistveno boljše predstave o sistemu ali bistveno natančnejšega odziva. Poveča pa zahtevnost računskih postopkov in čas reševanja.

Prostostne stopnje so med seboj linearno neodvisni nekonstitutivni parametri sistema, s katerimi je obnašanje sistema enolično določeno. Imenujemo jih tudi osnovne spremenljivke.

Število prostostnih stopenj je odvisno od:

• vrste osnovne fizikalne spremenljivke (skalar, vektor, tenzor),
• konstitucijskih lastnosti (sistem masnih točk, togo telo, deformabilno telo),
• prostorske razsežnosti (1D, 2D, 3D).

Število prostostnih stopenj materialne točke:

Spremenljivka	3D	2D	1D
skalar	1	1	1
vektor	3	2	1
tenzor 2. reda	9	4	1

Število prostostnih stopenj opredeljuje naravo sistema: diskretni sistemi imajo končno, kontinualni pa neskončno število prostostnih stopenj.

Za enoprostostni diskretni sistem je značilno:

• vselej ga je mogoče obravnavati v 1D prostoru,
• ne glede na značaj osnovne fizikalne spremenljivke ga je mogoče obravnavati skalarno z eno prostostno stopnjo,
• reševanje se prevede na reševanje ene same enačbe.

Za sistem zapišemo eno ravnotežno (gibalno) enačbo z eno neznanko, npr. $F = ku$ ali $m\frac{d^2x}{dt^2} = \sum_i F_i(x)$. Stacionarni problem da algebrajsko enačbo (AE), nestacionarni pa navadno diferencialno enačbo (NDE). Enačbo rešimo analitično ali numerično.

Za večprostostni diskretni sistem je značilno:

• spremenljivke, ki določajo prostost sistema, so med seboj linearno neodvisne,
• reševanje se prevede na reševanje sistema enačb: sistema algebrajskih enačb (SAE) pri stacionarnih in sistema NDE (SNDE) pri nestacionarnih problemih.

Ko število prostostnih stopenj preseže vse meje, preide diskretni sistem pri posebnih pogojih v kontinualni sistem.

Odvisna je od vrste osnovne fizikalne spremenljivke (skalar, vektor, tenzor), prostorske razsežnosti in časovne odvisnosti obravnavanega fizikalnega sistema.

Sistem	$\frac{d}{dt}() = 0$	$\frac{d}{dt}() \neq 0$
diskretni, enoprostostni	AE	NDE
diskretni, večprostostni	SAE	SNDE
kontinualni (skalar), 1D	NDE	PDE
kontinualni (skalar), 2D, 3D	PDE	PDE
kontinualni (vektor, tenzor)	SPDE	SPDE

AE – algebrajska enačba, NDE – navadna diferencialna enačba, PDE – parcialna diferencialna enačba, S – sistem.

Problem mora biti časovno neodvisen. Element mora biti raven, iz linearno elastičnega gradiva in obremenjen le s točkovno ali zvezno porazdeljeno osno obtežbo ter s temperaturno spremembo, ki je konstantna po prerezu.

Razmislek: Pri prečni obremenitvi napetost po prerezu ni več konstantna, zato problem ni več osni.

V vodilni enačbi je zajeto:

• statično ravnotežje vseh obremenitev na diferencialnem elementu Iz tega sledi $\frac{dN}{dx} = -n(x)$.
• deformacijska konsistentnost: deformacija je konsistentna s pomikom, $\frac{du}{dx} = \varepsilon_{xx}$, kjer je $\varepsilon_{xx} = \varepsilon_{xx}^{\sigma} + \varepsilon_{xx}^{T}$ (prispevek napetosti in temperature).
• konstitucijsko obnašanje: Hookov zakon, $\varepsilon_{xx}^{\sigma} = \frac{\sigma_{xx}}{E}$ in $\varepsilon_{xx}^{T} = \alpha \Delta T$.

Fizikalni spremenljivki sta:

• vzdolžni pomik $u(x)$ – primarna spremenljivka,
• notranja osna sila $N(x)$ – sekundarna spremenljivka, izražena s primarno: $N(x) = E A \left(\frac{du}{dx} - \alpha \Delta T\right)$.

Glede na robne pogoje sta konjugirani veličini.

Statično ravnotežje diferencialnega elementa $dx$: $$dN = -n\,dx \quad \Rightarrow \quad \frac{dN}{dx} = -n(x)$$ Deformacijska konsistentnost: $$\frac{du}{dx} = \varepsilon_{xx}, \quad \varepsilon_{xx} = \varepsilon_{xx}^{\sigma} + \varepsilon_{xx}^{T} \quad \Rightarrow \quad \varepsilon_{xx}^{\sigma} = \frac{du}{dx} - \varepsilon_{xx}^{T}$$ Konstitucijsko obnašanje (Hookov zakon): $$\varepsilon_{xx}^{\sigma} = \frac{\sigma_{xx}}{E}, \quad \varepsilon_{xx}^{T} = \alpha\Delta T, \quad \sigma_{xx} = \frac{N}{A} \quad \Rightarrow \quad N = E A \varepsilon_{xx}^{\sigma}$$ Iz tega sledi notranja osna sila: $$N = EA\left(\frac{du}{dx} - \alpha\Delta T\right)$$ Note: $E$, $A$, $\alpha$ in $\Delta T$ so lahko odvisni od $x$.

Ko to vstavimo v $\frac{dN}{dx} = -n(x)$, dobimo vodilno enačbo problema: $$\frac{d}{dx}\left[EA\left(\frac{du}{dx} - \alpha\Delta T\right)\right] = -n(x) \ ; \quad x \in [0, L]$$ To je navadna diferencialna enačba 2. reda z vzdolžnim pomikom $u(x)$ kot osnovno spremenljivko. V celoti opredeljuje spreminjanje $u(x)$ in $N(x)$ na 1D elementu.

Vodilna enačba upošteva le zvezno porazdeljeno obremenitev v polju elementa in temperaturno spremembo vzdolž elementa. Da je rešitev konsistentna tudi s pomiki in obremenitvami v krajiščih, mora zadostiti robnim pogojem na obeh krajiščih.

Robni pogoji so podani z znano vrednostjo primarne ali sekundarne spremenljivke v vsakem krajišču: $$x = 0: \quad u(0) = u_J \quad ali \quad N(0) = \left[EA\left(\frac{du}{dx} - \alpha\Delta T\right)\right]_{x=0} = -F_x^J$$ $$x = L: \quad u(L) = u_K \quad ali \quad N(L) = \left[EA\left(\frac{du}{dx} - \alpha\Delta T\right)\right]_{x=L} = +F_x^K$$ $u(x)$ in $N(x)$ sta konjugirani veličini. Na istem robu je ena izmed njiju znana, druga pa neznana.

Če so vsi parametri vodilne enačbe ($n(x)$, $\Delta T(x)$, $A(x)$, $E(x)$, $\alpha(x)$) na celotnem intervalu $x \in [0, L]$ vsak podani z enim funkcijskim predpisom, je tudi rešitev $u(x)$ podana z enim funkcijskim predpisom. Takrat so $u(x)$, $\frac{du}{dx}$ in $N(x)$ zvezne in zvezno odvedljive funkcije (1 polje).

Če ima vsaj en parameter na intervalu vsaj dva funkcijska predpisa, interval razdelimo na podintervale $[a_k, b_k]$ (več polj). Na meji med podintervaloma morajo biti izpolnjeni pogoji konsistentnosti prehoda. Ti so vedno fizikalno pogojeni, zato njihovo nespoštovanje vedno vodi do napačne rešitve.

Pogoja konsistentnosti prehoda sta:

• zveznost primarne spremenljivke (nerazdružljivost snovnih točk): $u_{k}(b_{k}) = u_{k+1}(a_{k+1})$
• statično ravnotežje notranjih sil in zunanje obremenitve na meji: $N_{k}(b_{k}) = N_{k+1}(a_{k+1}) + F_{k,k+1}$

$N(x)$ je na meji nezvezna le, če na meji deluje koncentrirana sila $F_{k,k+1}$. Drugi pogoj zapišemo s primarno spremenljivko: $$\left(EA\frac{du}{dx}\right)_k\bigg|_{x=b_k} - \left(EA\frac{du}{dx}\right)_{k+1}\bigg|_{x=a_{k+1}} = F_{k,k+1} + (EA\alpha\Delta T)_k\big|_{x=b_k} - (EA\alpha\Delta T)_{k+1}\big|_{x=a_{k+1}}$$

Eksaktni funkcijski predpis $u(x)$ dobimo analitično z dvojnim nedoločenim integriranjem vodilne enačbe. Ker je enačba NDE 2. reda, dobimo 2 integracijski konstanti, ki ju določimo iz robnih pogojev.

Pri več podintervalih dobimo 2 konstanti na podinterval. Določimo jih iz robnih pogojev in pogojev konsistentnosti prehoda.

Rešitev je eksaktna, ko zadošča vodilni enačbi v vsaki točki območja $x \in [0, L]$ ter robnim pogojem (in pogojem konsistentnosti prehoda).

Rešitev je aproksimativna, ko zadošča robnim pogojem in pogojem konsistentnosti prehoda, vodilni diferencialni enačbi pa ne v vsaki točki območja $x \in [0, L]$, temveč le v izbranih točkah oz. podobmočjih.

Aproksimativna rešitev je zasnovana na končni množici parametrov $c_i$: $$u_N(x) = f(x, c_i \ ; \ i = 0, 1, \dots, N), \qquad \lim_{N \to \infty} u_N(x) = u(x)$$

Pri funkcijskem pristopu je aproksimativna rešitev zasnovana na končni množici $\{ \Psi_{i} (x) \}$ izbranih aproksimacijskih funkcij iste družine. Funkcije $\Psi_i(x)$ so zvezne in zvezno odvedljive.

Aproksimativna rešitev $u_{N}(x)$ je: $$u_{N}(x) = \sum_{i = 0}^{N}c_{i} \Psi_{i} (x)\ ;\ x\in [0, L]$$ V funkcijskem smislu jo opredeljuje izbira funkcij $\Psi_i(x)$, po vrednosti pa koeficienti $c_i$. Neznane koeficiente $c_i$ določimo v postopku reševanja.

Vsako zvezno in poljubnokrat odvedljivo funkcijo lahko razvijemo v Taylorjevo (potenčno) vrsto okoli $x_0$. Zato je smiselna izbira potenc: $$\Psi_i(x) = (x - x_0)^i \ ; \quad i = 0, 1, 2, \dots, N$$ Da nabor $\Psi_{i}(x)$ omogoča popis eksaktne rešitve $u(x)$, mora biti kompleten. Vsebovati mora po vrsti vse potence od najnižje ($i = 0$) do najvišje ($i = N$). Izpustitev ene izmed nižjih potenc ni popravljiva.

$u_N(x)$ ima $N+1$ neznanih koeficientov $c_i$, zato potrebujemo sistem $N+1$ linearno neodvisnih enačb. Da je aproksimativna rešitev fizikalno konsistentna in verodostojna, mora zadostiti ključnim enačbam problema:

• enačbam robnih pogojev,
• enačbam pogojev konsistentnosti prehoda,
• in v čim večji meri vodilni diferencialni enačbi. Iz nje dobimo manjkajoče enačbe, ko jo zapišemo v izbranih točkah.

Pri več podintervalih $u_N^{(k)}(x)$ iščemo za vsak podinterval posebej. Takrat potrebujemo $\sum_{k=1}^{n}(N_k + 1)$ enačb.

Pri interpolacijskem pristopu je aproksimativna rešitev zasnovana na končni množici $\{ c_i\}$ diskretnih parametrov. Ti aproksimirajo vrednosti primarne spremenljivke v izbranih točkah območja: $c_i \equiv u_i \approx u(x_i)$.

Aproksimativno rešitev zapišemo kot: $$u_N(x) = \Psi_{int}(x \times \{c_i\ ;\ i = 0,1,\dots,N\})\ ;\ x\in[0,L]$$ V funkcijskem smislu je opredeljena šele z izbiro interpolacijske funkcije $\Psi_{int}$. Ta ekstrapolira diskretne vrednosti $u_i$ na celotno območje $[0, L]$. Pri istem naboru $\{u_i\}$ dobimo z različnimi interpolacijskimi funkcijami različne aproksimacije.

Najprej moramo torej določiti neznane vrednosti $u_i$ v točkah $x_i$. Za $N+1$ neznank potrebujemo sistem $N+1$ linearno neodvisnih enačb. Da je rešitev fizikalno smiselna, sistem tvorimo po vrsti iz robnih pogojev, pogojev konsistentnega prehoda in vodilne DE.

Aproksimacijo gradimo na diskretnih vrednostih $u_i$. Fizikalna konsistenca pa je pogojena z diferencialnimi zvezami (vodilna enačba, $N(x)$), ki zahtevajo funkcijsko obravnavo. Zato diferencialne zveze pretvorimo v diferenčne.

Če je $u(x)$ v okolici $x_0$ zvezna in zvezno odvedljiva, jo razvijemo v Taylorjevo vrsto: $$u(x_0 + h) = u(x_0) + \frac{h^1}{1!}\frac{d^1u}{dx^1}(x_0) + \frac{h^2}{2!}\frac{d^2u}{dx^2}(x_0) + \frac{h^3}{3!}\frac{d^3u}{dx^3}(x_0) + \ ...$$ Tako izrazimo vrednost v sosednjih točkah $x_0 + h^+$ in $x_0 - h^-$:

Vpeljemo oznake: $$u_0 = u(x_0) \ ;\ u^+ = u(x_0 + h^+)\ ;\ u^- = u(x_0 - h^-)$$ $$u^+ = u_0 + \frac{(h^+)^1}{1!}\frac{d^1u_0}{dx^1} + \frac{(h^+)^2}{2!}\frac{d^2u_0}{dx^2} + \frac{(h^+)^3}{3!}\frac{d^3u_0}{dx^3} + \ ...$$ $$u^- = u_0 - \frac{(h^-)^1}{1!}\frac{d^1u_0}{dx^1} + \frac{(h^-)^2}{2!}\frac{d^2u_0}{dx^2} - \frac{(h^-)^3}{3!}\frac{d^3u_0}{dx^3} + \ ...$$ Če velja $h^3 \ll h < 1$, zanemarimo člene s 3. in višjimi odvodi.

Z eliminacijo 2. odvoda dobimo prvi, z eliminacijo 1. odvoda pa drugi odvod: $$\frac{du_0}{dx} \approx \frac{1}{h^+ + h^-}\left[\frac{h^-}{h^+}(u^+ - u_0) - \frac{h^+}{h^-}(u^- - u_0)\right]$$ $$\frac{d^2u_0}{dx^2} \approx \frac{2}{h^+ + h^-}\left[\frac{1}{h^+}(u^+ - u_0) + \frac{1}{h^-}(u^- - u_0)\right]$$ Odvode tako nadomestimo z diferenčnimi operatorji: $\frac{d^r u_0}{dx^r} \approx D^r u_0$.

Centralno shemo dobimo, ko sta točki enako oddaljeni: $h^- = h^+ = h$.

Za 1. odvod $u^+$ in $u^-$ odštejemo, da se eliminira 2. odvod: $$u^+ - u^- \approx \left(u_0 + h\frac{du_0}{dx} + \frac{h^2}{2}\frac{d^2u_0}{dx^2}\right) - \left(u_0 - h\frac{du_0}{dx} + \frac{h^2}{2}\frac{d^2u_0}{dx^2}\right) = 2h\frac{du_0}{dx}$$ Iz tega sledi: $$\frac{du_0}{dx} \approx \frac{u^+ - u^-}{2h} = D^1 u_0$$

Za 2. odvod $u^+$ in $u^-$ seštejemo, da se eliminira 1. odvod: $$u^+ + u^- \approx \left(u_0 + h\frac{du_0}{dx} + \frac{h^2}{2}\frac{d^2u_0}{dx^2}\right) + \left(u_0 - h\frac{du_0}{dx} + \frac{h^2}{2}\frac{d^2u_0}{dx^2}\right) = 2u_0 + h^2\frac{d^2u_0}{dx^2}$$ Iz tega sledi: $$\frac{d^2u_0}{dx^2} \approx\frac{u^+ - 2u_0 + u^-}{h^2} = D^2 u_0$$ (note: $u^-$, $u_0$ in $u^+$ za $i$-to točko pišemo tudi kot $u_{i-1}$, $u_i$ in $u_{i+1}$.)

Metoda končnih razlik (MKR) je aproksimativna metoda, ki temelji na interpolacijskem pristopu. Sistem enačb za neznane vrednosti $u_k$ temelji na pretvorbi diferencialnih operatorjev v diferenčne (centralne razlike $D^1$ in $D^2$).

Postopek:

1. Območje razdelimo na $N$ podintervalov s korakom $h$. Neznanke so vrednosti $u_k \approx u(x_k)$ v $N+1$ točkah.
2. Sistem linearno neodvisnih enačb tvorimo tako, da v čim večji meri izpolnimo:
   • robne pogoje,
   • pogoje konsistentnosti prehoda,
   • območno enačbo problema (vodilno DE v diferenčni obliki, npr. $EA\frac{u_{k+1} - 2u_k + u_{k-1}}{h^2} = -n_0$ v notranjih točkah).
3. Rešimo sistem. Iz $u_k$ izračunamo še $N_k$ z $D^1$.

Več točk (manjši $h$) načeloma da natančnejšo rešitev.

Točke vedno izberemo tudi na robu območja. Iz robnih pogojev tvorimo največ toliko enačb, kolikor je robnih pogojev.

• Robni pogoj s primarno spremenljivko zapišemo neposredno, npr. $u(0) = 0 \Rightarrow u_0 = 0$.
• Robni pogoj s sekundarno spremenljivko ($N$) vsebuje odvod, ki ga nadomestimo z $D^1$. Centralna razlika v robni točki potrebuje dodatno točko izven območja, npr. $D^1 u_4 = \frac{u_A - u_3}{2h}$.

Dodatna točka poveča število neznank za ena. Manjkajočo enačbo dobimo z zapisom vodilne DE tudi v robni točki ($D^2 u_4$). Vrednost v dodatni točki nima fizikalnega pomena.

Prvi pogoj $u_1(L_1) = u_2(L_1)$ je samodejno izpolnjen, ker točko postavimo na mejo med podintervaloma. Tam je ena sama neznana diskretna vrednost.

Drugi pogoj $N_1(L_1) - N_2(L_1) = 0$ vsebuje odvoda: $$EA_1\left(\frac{du_1}{dx}\right)\bigg|_{x=L_1} - E A_2\left(\frac{du_2}{dx}\right)\bigg|_{x = L_1} = EA_1\alpha\Delta T - EA_2\alpha\Delta T$$ $D^1$ sme za vsak podinterval uporabiti le vrednosti iz tega podintervala. Zato za vsak podinterval posebej dodamo dodatno točko ($B$ za 1. in $C$ za 2. podinterval). V diskretni obliki je zapis: $$A_1\frac{u_B-u_2}{2h_1} - A_2\frac{u_4 - u_C}{2h_2} = \alpha\Delta T(A_1 - A_2)$$ Za dodatni neznanki $u_B$ in $u_C$ zapišemo vodilno DE v mejni točki za vsak podinterval posebej.

(note: če je eksaktna rešitev polinom največ 2. stopnje, sta $D^1$ in $D^2$ eksaktna. Takrat je MKR rešitev v diskretnih točkah eksaktna.)

Reševanje vodilne diff. enačbe lahko prevedemo na reševanje ustrezne integralske enačbe. Izhodišče predstavlja že poznana vodilna enačba problema : $$\frac{d}{dx}[EA(\frac{du}{dx} - \alpha\Delta T)] = -n\space;\space x\in[0, L]$$ Simbolno lahko enačbo zapišemo kot $$D^2u(x) = f(x)$$ kjer je $$D^2 = \frac{d}{dx}[EA(\frac{d}{dx} - \alpha\Delta T)]\space\text{in}\space f(x) = -n(x)$$ Vodilno DE preoblikujemo v obliko: $$D^2u(x) - f(x) = 0$$ Enačbo nato pomnožimo s poljubno na območju $x\in[0, L]$ odvedljivo funkcijo $v(x)$ :$$[D^2u(x) - f(x)]\space v(x) = 0$$ Ker je izraz v oklepaju na celotnem območju $x\in[0,L]$ ničen, je ničen tudi produkt, ne glede na to kakšna je funkcija $v(x)$. Produkt integriramo po celotnem območju: $$\int_{0}^{L}[D^2u(x) - f(x)]\space v(x)dx = 0$$ Integrand je povsod ničen, zato je ničen tudi integral. Tako dobimo osnovno obliko integralske formulacije.

Za osno obremenjeni konstrukcijski element lahko to integralsko enačbo zapišemo kot : $$\int_{0}^{L}\frac{d}{dx}[E A(\frac{du}{dx} - \alpha\Delta T)]\space v(x) dx + \int_{0}^{L}n(x)v(x)dx = 0$$

Integralska formulacija zajema vse točke obravnavanega območja (to ne pomeni, da je rešitev eksaktna!). Pri šibki in inverzni obliki je red diferencialnega operatorja nad $u(x)$ nižji, kar zmanjša zahteve glede aproksimacijske funkcije.

Integralska formulacija je osnova za:

• Metodo končnih razlik (osnovna oblika)
• Metodo končnih elementov (šibka oblika)
• Metodo robnih elementov (inverzna oblika)

Integralov ni mogoče izvrednotiti brez funkcijske aproksimacije primarne spremenljivke: $$\hat{u}(x) = \sum_{k = 0}^N a_k \hat{\Psi}_k(x)$$ Neznanke so diskretne vrednosti $\{a_k\}$ (in robne vrednosti sekundarne spremenljivke). Robni pogoji dajo le del enačb (pri DE 2. reda dve).

Manjkajoče enačbe dobimo z izbiro poljubne funkcije $v(x)$. Vsaka izbrana funkcija $v_k(x)$ da eno enačbo:

• MKR: $v_k(x) = \delta(x-x_k)$ (Diracova funkcija v posameznih točkah),
• MKE: $v_k(x) = \hat{\Psi}_k(x)$ (Galerkinov pristop, iste funkcije kot pri aproksimaciji),
• MRE: $v_k(x)$ je rešitev enačbe $\tilde{D}^2v_k(x) = \delta(x-x_k)$.

Primer, 2-vozliščni KE: neznanke so $U_1^e, U_2^e, N_1^e, N_2^e$. Dve enačbi dobimo iz robnih pogojev v vozliščih, dve pa z izbiro $v_1 = \hat{\Psi}_0^e = 1-\frac{x_e}{L_e}$ in $v_2 = \hat{\Psi}_1^e = \frac{x_e}{L_e}$.

Izpeljava se začne pri osnovni obliki integralske formulacije, ki smo jo že izpeljali : $$\int_0^L[D^2u(x) - f(x)]\space v(x)\space dx = 0$$ Enačbo lahko preuredimo v naslednjo obliko : $$\int_0^LD^2u(x)v(x)dx = \int_0^Lf(x)v(x)dx$$ Če privzamemo, da je $v(x)$ vsaj enkrat odvedljiva funkcija, lahko integral na levi strani enačbe enkrat integriramo z uporabo Per-Partes metode $\int_{x_1}^{x_2}u(x)dv = uv|_{x1}^{x_2} - \int_{x_1}^{x_2}v(x)du$ : $$\int_0^LD^2 u(x)v(x)dx = D^1u(x)v(x)\biggr{|}_0^L - \int_0^LD^1u(x)\frac{dv(x)}{dx}dx$$ Red diferencialnega operatorja nad primarno spremenljivko $u(x)$ se je zmanjšal za 1: $D^2u(x) \rightarrow D^1u(x)$.

Celotno enačbo lahko sedaj zapišemo kot : $$\int_0^LD^1u(x)\frac{dv(x)}{dx}dx = D^1u(x)v(x)\biggr{|}_0^L - \int_0^Lf(x)v(x)dx$$ Tak zapis integralske enačbe imenujemo ŠIBKA OBLIKA INTEGRALSKE FORMULACIJE. Funkcija $v(x)$ mora izpolnjevati $\frac{dv}{dx} \neq 0$.

Za enoosno obremenjeni konstrukcijski element je diff. operator $D^1$ enak : $$D^1 = E(x)A(x)(\frac{d}{dx} - \alpha(x)\Delta T(x))$$ To lahko vstavimo v zgornjo enačbo in zapišemo : $$\int_0^L E A\frac{du}{dx}\frac{dv}{dx}dx = [E A(\frac{du}{dx} - \alpha\Delta T)]v(x)\biggr{|}_0^L - \int_0^Lf(x)v(x)dx + \int_0^LE A \alpha\Delta T \frac{dv}{dx}dx$$

Robne vrednosti sekundarne spremenljivke $N(0)$ in $N(L)$. Ob upoštevanju $N = EA(\frac{du}{dx} - \alpha\Delta T)$ se šibka oblika glasi: $$\int_0^LEA\frac{du}{dx}\frac{dv}{dx}dx = N(L)v(L) - N(0)v(0) - \int_0^Lf(x)v(x)dx + \int_0^LE A \alpha\Delta T \frac{dv}{dx}dx$$ Če je na robu podan robni pogoj za sekundarno spremenljivko, sta $N(0)$ oziroma $N(L)$ znani vrednosti. Robnih vrednosti primarne spremenljivke $u(0)$, $u(L)$ šibka oblika ne vsebuje.

DE je izpolnjena le v integralskem smislu po celotnem območju, ne v vsaki točki. Z aproksimativno rešitvijo $\hat{u}(x)$ DE ni izpolnjena v vseh točkah območja, zato rešitev ni eksaktna. Pri MKR ($v_k = \delta(x-x_k)$) je DE izpolnjena le v izbranih točkah $x_k$.

Zapišimo šibko obliko integralske formulacije za 1D osno obremenjeni element : $$\int_0^L\tilde{D}^1u(x)\frac{dv(x)}{dx}dx = N(x)v(x)\biggr{|}_0^L - \int_0^Lf(x)v(x)dx + \int_0^Lf_T(x) \frac{dv}{dx}dx$$ kjer je $\tilde{D}^1 = E(x)A(x)(\frac{d}{dx})$ in je $f_T(x) = E(x)A(x)\alpha(x)\Delta T(x)$.

S takim zapisom lahko vidimo, da lahko levo stran integralske formulacije še enkrat Per-Partes integriramo . $$\int_0^L \tilde{D}^1 u(x) \frac{dv(x)}{dx}dx = u(x) \tilde{D}^1v(x)\biggr{|}_0^L - \int_0^Lu(x)\tilde{D}^2v(x) dx$$ S tem integralska formulacija dobi sledečo obliko : $$\int_0^Lu(x)\tilde{D}^2v(x)dx =u(x)\tilde{D}^1v(x)|_0^L - N(x)v(x)\biggr{|}_0^L + \int_0^Lf(x)v(x)dx - \int_0^Lf_T(x)\frac{dv(x)}{dx}dx$$ Dobljeni obliki integralske enačbe rečemo INVERZNA OBLIKA INTEGRALSKE FORMULACIJE. Diferencialni operator je v celoti prešel s primarne spremenljivke na funkcijo $v(x)$ ($D^2u \rightarrow \tilde{D}^2v$, kjer je $\tilde{D}^2 = \frac{d}{dx}(EA\frac{d}{dx})$), zato mora veljati $\frac{d^2v}{dx^2} \neq 0$.

Inverzna oblika vsebuje robne vrednosti primarne spremenljivke $u(0)$ in $u(L)$ ter sekundarne spremenljivke $N(0)$ in $N(L)$. Pri podanem robnem pogoju so to znane vrednosti.

Vozlišča KE so podana s koordinatami v globalnem K.S. $X$. Te koordinate določajo lego in velikost KE (npr. $L_1 = X_2 - X_1$). Vsak KE ima še svoj lokalni K.S. $x_e\in[0, L_e]$, v katerem zapišemo aproksimacijo.

Aproksimacijo izvedemo ločeno na vsakem podobmočju $x_e\in[0,L_e]$, imenovanem končni element (KE), ki ima $N_e$ vozlišč: $$\hat{u}_e = \hat{u}_e(x_e) = \sum_{k=0}^{N_e-1}a_k^e\hat{\Psi}_k^e(x_e) \space ; \space x_e\in[0,L_e]$$ Interpolacijske funkcije $\hat{\Psi}_k^e$ so Lagrangevi polinomi, ki izpolnjujejo pogoj (Kroneckerjev delta): $$\hat{\Psi}_k^e(x_i) = \delta_{ik} = \begin{cases}1, &\text{za } x_i = x_k\\0, &\text{za }x_i \neq x_k\end{cases}$$ Zato so koeficienti kar vrednosti primarne spremenljivke v vozliščih: $a_k^e = \hat{u}_e(x_k)$. Število vozlišč $N_e$ določa stopnjo polinoma ($N_e - 1$).

Aproksimacijo vstavimo v šibko obliko integralske formulacije : $$\int_0^{L_e} D^1\hat{u}(x) \frac{dv(x)}{dx}dx = N(x)v(x)\biggr{|}_0^{L_e} - \int_0^{L_e}f(x)v(x)dx + \int_0^{L_e}f_T(x)\frac{dv(x)}{dx}dx$$ Funkcije $v(x)$ izberemo po Galerkinu: $v_k = \hat{\Psi}_k^e$. Za 2-vozliščni KE glej naslednje vprašanje.

V primeru 2-vozliščnega KE je aproksimacija primarne spremenljivke $\hat{u}_e(x_e)$ po njegovem območju zasnovana na 2 polinomih prvega reda:

S tem je aproksimativna rešitev linearna funkcija, določena z 2 konstantama : $$\hat{u}_e(x_e)= a_0^e(1- \frac{x_e}{L_e}) + a_1^e(\frac{x_e}{L_e})$$

Pri $x_e = 0$ in $x_e = L_e$ velja : $$\hat{u}_e(0) = a_0^e = U_1^e\space,\space\hat{u}_e(L_e) = a_1^e = U_2^e$$ Konstanti sta torej pomika v vozliščih KE

Predpostavimo, da je KE:

• konstantnega prereza ($A(x_e) = A_0$)
• konstantnih materialnih lastnosti ($E(x_e) = E_0$ in $\alpha(x_e) = \alpha_0$)

Šibka oblika integralske formulacije : $$E_0A_0\int_0^{L_e}\frac{d\hat{u}_e}{dx_e}\frac{dv}{dx_e}dx_e = N_e(L_e)v(L_e) - N_e(0)v(0) + \int_0^{L_e}n(x_e)v(x_e)dx_e + \int_0^{L_e} E_0A_0\alpha_0\Delta T \frac{dv}{dx_e}dx_e$$ $N(0)$ in $N(L_e)$ sta vozliščni vrednosti notranje osne sile: $$N_e(0) = N_1^e\space;\space N_e(L_e) = N_2^e$$

Vstavimo aproksimacijo, za katero je $\frac{d\hat{u}_e}{dx_e} = U_1^e(-\frac{1}{L_e}) + U_2^e(\frac{1}{L_e})$ : $$\begin{aligned}&E_0A_0[U_1^e(-\frac{1}{L_e}) + U_2^e(\frac{1}{L_e})]\int_0^{L_e}\frac{dv}{dx_e}dx_e = \\ &= N_2^ev(L_e) - N_1^ev(0) + \int_0^{L_e}n(x_e)v(x_e)dx_e + \int_0^{L_e} E_0A_0\alpha_0\Delta T \frac{dv}{dx_e}dx_e\end{aligned}$$

V enačbi so 4 neznanke: $U_1^e$, $U_2^e$, $N_1^e$, $N_2^e$. 2 enačbi (vodilna enačba je DE 2. reda) izhajata iz poznanih vrednosti primarne ali sekundarne spremenljivke v vozliščih KE.

Manjkajoči 2 enačbi dobimo z izbiro poljubne funkcije $v(x_e)$. Po Galerkinovem pristopu izberemo funkciji, uporabljeni v aproksimaciji $\hat{u}_e(x_e)$ : $$v_1(x_e) = \hat{\Psi}_0^e(x_e) = 1 - \frac{x_e}{L_e}$$ $$v_2(x_e) = \hat{\Psi}_1^e(x_e) = \frac{x_e}{L_e}$$

1. enačba ($v = v_1$, $\frac{dv_1}{dx_e} = -\frac{1}{L_e}$, $v_1(0) = 1$, $v_1(L_e) = 0$) : $$\frac{E_0A_0}{L_e}[U_1^e - U_2^e] = -N_1^e + \int_0^{L_e}n(1-\frac{x_e}{L_e})dx_e + \int_0^{L_e}E_0A_0\alpha_0\Delta T(-\frac{1}{L_e})dx_e$$
2. enačba ($v = v_2$, $\frac{dv_2}{dx_e} = \frac{1}{L_e}$, $v_2(0) = 0$, $v_2(L_e) = 1$) : $$\frac{E_0A_0}{L_e}[-U_1^e + U_2^e] = N_2^e + \int_0^{L_e}n(\frac{x_e}{L_e})dx_e + \int_0^{L_e}E_0A_0\alpha_0\Delta T(\frac{1}{L_e})dx_e$$ Matrična oblika : $$\frac{E_0A_0}{L_e}\begin{bmatrix} 1 &-1 \\ -1 &1 \end{bmatrix} \begin{pmatrix} U_1^e \\ U_2^e \end{pmatrix} = \begin{pmatrix} -N_1^e \\ N_2^e \end{pmatrix} + \begin{pmatrix} F_1^e \\ F_2^e \end{pmatrix}$$ kjer je $F_i^e = F_{in}^e + F_{iT}^e$ (vprašanji 12 in 13). Togostna matrika KE je simetrična in odvisna le od geometrije in materiala KE.

Kot ekvivalentni sili v vozliščih KE : $$F_{1n}^e = \int_0^{L_e}n(1-\frac{x_e}{L_e})dx_e$$ in $$F_{2n}^e = \int_0^{L_e}n(\frac{x_e}{L_e})dx_e$$ Za konstantno $n = n_0$ je $F_{1n}^e = F_{2n}^e = \frac{n_0L_e}{2}$: rezultanta $n_0L_e$ se razdeli na obe vozlišči.

Kot ekvivalentni sili v vozliščih, ki povzročita enak raztezek KE kot sprememba temperature : $$F_{1T}^e = \int_0^{L_e}E_0A_0\alpha_0\Delta T(x_e)(-\frac{1}{L_e})dx_e$$ in $$F_{2T}^e = \int_0^{L_e}E_0A_0\alpha_0\Delta T(x_e)(\frac{1}{L_e})dx_e$$ Za konstantno $\Delta T = \Delta T_0$ je $F_{1T}^e = -E_0A_0\alpha_0\Delta T_0$ in $F_{2T}^e = +E_0A_0\alpha_0\Delta T_0$.

1. Geometrijski model: območje razdelimo na KE, vozlišča podamo s koordinatami v globalnem K.S.
2. Na vsakem KE aproksimiramo primarno spremenljivko $\hat{u}_e(x_e) = \sum_{k=0}^{N_e-1} a_k^{e} \hat{\Psi}_k^{e} (x_e)$ (koeficienti $a_k^e$ so vozliščni pomiki) in jo vstavimo v šibko obliko integralske formulacije.
3. Funkcije $v(x_e)$ izberemo po Galerkinu ($v_k = \hat{\Psi}_k^e$). Za 2-vozliščni KE dobimo sistem KE: $$\frac{E_0A_0}{L_e}\begin{bmatrix} 1 &-1 \\ -1 &1 \end{bmatrix} \begin{Bmatrix} U_1^e \\ U_2^e \end{Bmatrix} = \begin{Bmatrix} -N_1^e \\ N_2^e \end{Bmatrix} + \begin{Bmatrix} F_1^e \\ F_2^e \end{Bmatrix}$$
4. Izračunamo obremenitve KE $F_i^e = F_{in}^e + F_{iT}^e$ (porazdeljena in temperaturna obremenitev).
5. Sisteme KE sestavimo v globalni sistem (pri več KE ob upoštevanju konsistentnosti prehoda).
6. Upoštevamo robne pogoje: v vsakem robnem vozlišču je znan pomik ali sila.
7. Rešimo sistem in dobimo neznane pomike in sile v vozliščih.
8. Po območju KE izračunamo $\hat{u}_e(x_e)$, deformacijo $\varepsilon_{xx} = \frac{d\hat{u}_e}{dx_e}$, osno silo $N = E_0A_0(\frac{d\hat{u}_e}{dx_e} - \alpha_0\Delta T)$ in napetost $\sigma_{xx} = E_0(\frac{d\hat{u}_e}{dx_e} - \alpha_0\Delta T)$. Pri 2-vozliščnem KE so te veličine po KE konstantne.

Pogoje konsistentnega prehoda upoštevamo za primarno in za sekundarno spremenljivko v vozlišču, ki je skupno dvema KE.

Primarna spremenljivka (zveznost pomika): skupno vozlišče ima en sam pomik. Ker globalna os $X$ in lokalni osi $x_e$ sovpadata, velja : $$U_1^1 = U_1\space;\space U_2^1 = U_2^2 = U_2\space;\space U_3^2 = U_3$$ Problem ima tako 3 neznane vrednosti primarne spremenljivke. Sistem vsakega KE razširimo na vse 3 globalne neznanke:

KE 1 : $$\frac{E_0A_0}{L_1}\begin{bmatrix} 1&-1&0\\-1&1&0\\0&0&0 \end{bmatrix}\begin{Bmatrix}U_1\\U_2\\U_3\end{Bmatrix} = \begin{Bmatrix} -N_1^1\\N_2^1\\0\end{Bmatrix} + \begin{Bmatrix} F_{1n}^1\\F_{2n}^1\\0\end{Bmatrix}$$ KE 2 : $$\frac{E_0A_0}{L_2}\begin{bmatrix}0&0&0\\0&1&-1\\0&-1&1\end{bmatrix}\begin{Bmatrix}U_1\\U_2\\U_3\end{Bmatrix} = \begin{Bmatrix} 0\\-N_2^2\\N_3^2\end{Bmatrix} + \begin{Bmatrix} 0\\F_{2n}^2\\F_{3n}^2\end{Bmatrix}$$ Ker v obeh sistemih nastopajo iste neznanke, ju seštejemo : $$E_0A_0\begin{bmatrix}\frac{1}{L_1}&-\frac{1}{L_1}&0\\-\frac{1}{L_1}&\frac{1}{L_1}+\frac{1}{L_2}&-\frac{1}{L_2}\\0&-\frac{1}{L_2}&\frac{1}{L_2}\end{bmatrix}\begin{Bmatrix}U_1\\U_2\\U_3\end{Bmatrix} = \begin{Bmatrix} -N_1^1\\N_2^1-N_2^2\\N_3^2\end{Bmatrix} + \begin{Bmatrix}F_{1n}^1\\F_{2n}^1 + F_{2n}^2\\F_{3n}^2\end{Bmatrix}$$ Ta sistem 3 enačb se nanaša na celotno območje in ima 7 neznank ($U_1,U_2,U_3, N_1^1, N_2^1, N_2^2, N_3^2$).

Robna pogoja $U_1 = 0$ in $N_3^2 = 0$ odpravita 2 neznanki.

Sekundarna spremenljivka (ravnotežje sil v vozlišču 2, $dx_1, dx_2 \rightarrow 0$) : $$-N_2^1 + N_2^2 + F_0 - n_0dx_1 - n_0dx_2 = 0 \rightarrow N_2^1 - N_2^2 = F_0$$

Pogoj konsistentnega prehoda za sekundarno spremenljivko je torej statično ravnotežje notranjih in zunanjih sil na meji med KE. Neznani sili $N_2^1$ in $N_2^2$ smo nadomestili z znano silo $F_0$. Ostanejo 3 neznanke ($U_2,U_3,N_1^1$) in sistem je rešljiv (za $n(x) = n_0$): $$E_0A_0\begin{bmatrix}\frac{1}{L_1}&-\frac{1}{L_1}&0\\-\frac{1}{L_1}&\frac{1}{L_1}+\frac{1}{L_2}&-\frac{1}{L_2}\\0&-\frac{1}{L_2}&\frac{1}{L_2}\end{bmatrix}\begin{Bmatrix}0\\U_2\\U_3\end{Bmatrix} = \begin{Bmatrix} -N_1^1\\F_0\\0\end{Bmatrix} + \begin{Bmatrix}-\frac{n_0L_1}{2}\\-\frac{n_0(L_1+L_2)}{2}\\-\frac{n_0L_2}{2}\end{Bmatrix}$$