Metóda konečných prvkov je numerická metóda slúžiaca na simuláciu priebehov napätia, deformácii, vlastných frekvencií, prúdenie tepla, javov elektromagnetizmu, prúdenie tekutín a pod. na vytvorenom fyzikálnom modeli. Sme budeme riešiť diferenciálnu rovnicu 2. rádu v 1D v stacionárnom prípade. \[\begin{equation}
- \frac{d}{dx} \left( a(x) \frac{du}{dx} \right) = f(x), \quad x \in (0, L),\tag{1}
\label{strong_form}
\end{equation}\]
kde neznámou je funkcia \(u(x)\). Vstupmi (dátami) úlohy sú materiálová charakteristika \(a(x)\), intenzita zdrojov (resp. vonkajších síl) \(f(x)\) a okrajové podmienky zadané na hranici výpočtovej oblasti \(\Omega = (0, L)\) (čiže v bodoch \(x = 0\) a \(x = L\) ).
Funkcii $u(x) $ sa niekedy hovorí primárna neznáma a funkcii
kde \(\alpha_1, \beta_1, g_1\) a \(\alpha_2, \beta_2, g_2\) sú zadané konštanty.
Prípad \(\beta_i = 0\) zodpovedá Dirichletovej okrajovej podmienke, kde zadávame hodnotu primárnej neznámej \(u(x)\) na hranici. Ak \(\alpha_i = 0\), hovoríme o Neumannovej okrajovej podmienke, pri ktorej je vlastne na hranici zadaná hodnota sekundárnej neznámej \(q(x)\). Ak je \(\alpha_i\) aj \(\beta_i\) nenulové, ide o Newtonovu okrajovú podmienku.
2 Diskretizácia výpočtovej oblasti
Výpočtovú oblasť \(\Omega = (0, L)\) diskretizujeme na \(N\) elementov \(\Omega^e, e = 1, \dots, N\). Každý element je interval \(\Omega^e = [x_A^e, x_B^e]\), pričom intervaly môžu mať rôznu veľkosť \(h^e = |\Omega^e| = x_B^e - x_A^e.\)
Približné (numerické) riešenie na elemente \(\Omega^e\) budeme označovať ako \(U^e(x), x \in \Omega^e\). Celkové približné riešenie \(U(x), x \in [0, L]\) získame poskladaním funkcií \(U^e(x)\) do jednej funkcie, čiže
\[
U(x) =
\begin{cases}
U^1(x), & x \in \Omega^1, \\
U^2(x), & x \in \Omega^2, \\
\vdots & \\
U^N(x), & x \in \Omega^N.
\end{cases}
\]
3 Sústava rovníc na elemente
3.1 Slabá formulácia
Vyberme jeden element \(\Omega^e = [x_A^e, x_B^e]\). Vynásobme diferenciálnu rovnicu \(\eqref{strong_form}\) váhovou funkciou \(w(x)\) a výsledok zintegrujme na elemente \(\Omega^e\). Dostaneme variačnú formuláciu:
\[\begin{equation}\int_{x_A^e}^{x_B^e} -\frac{d}{dx} \left( a \frac{du}{dx} \right) w \, dx = \int_{x_A^e}^{x_B^e} f w \, dx. \tag{2}
\label{eq2}
\end{equation}\]
Na ľavej strane prenesieme jednu deriváciu na váhovú funkciu, využijeme na to integráciu per partes:
\[\int_a^b g' h \, dx = [g h]_a^b - \int_a^b g h' \, dx, \]
pričom \(g' = -\frac{d}{dx} \left( a \frac{du}{dx} \right)\), \(h = w\), a teda \(g = -a \frac{du}{dx}\), \(h' = \frac{dw}{dx}\).
Dostaneme:
\(\left[ -a \frac{du}{dx} w \right]_{x_A^e}^{x_B^e} + \int_{x_A^e}^{x_B^e} a \frac{du}{dx} \frac{dw}{dx} \, dx = \int_{x_A^e}^{x_B^e} f w \, dx.\)
A po rozpísaní hraničného člena máme: \[\begin{equation}
- \left( a \frac{du}{dx} \Big|_{x = x_B^e} \right) w(x_B^e) - \left( -a \frac{du}{dx} \Big|_{x = x_A^e} \right) w(x_A^e) + \int_{x_A^e}^{x_B^e} a \frac{du}{dx} \frac{dw}{dx} \, dx = \int_{x_A^e}^{x_B^e} f w \, dx.
\tag{3}
\label{eq3}
\end{equation}\]
Zaveďme označenia: \[
Q_A^e := \left( -a \frac{du}{dx} \Big|_{x = x_A^e} \right) \quad \text{tok dovnútra elementu zľava},
\]\[\begin{equation}
Q_B^e := \left( a \frac{du}{dx} \Big|_{x = x_B^e} \right) \quad \text{tok dovnútra elementu sprava}.
\tag{4}
\label{eq4}
\end{equation}\]
Hraničné členy v rovnici \(\eqref{eq3}\) prehodíme na pravú stranu a dostávame slabú formuláciu diferenciálnej rovnice \(\eqref{strong_form}\) na elemente \(\Omega^e\):
\[\begin{equation}
\int_{x_A^e}^{x_B^e} a \frac{du}{dx} \frac{dw}{dx} \, dx = \int_{x_A^e}^{x_B^e} f w \, dx + Q_A^e w(x_A^e) + Q_B^e w(x_B^e). \tag{5}
\label{eq5}
\end{equation}\] ## Aproximačné funkcie
Približné riešenie \(U^e(x)\) na elemente \(\Omega^e\) budeme hľadať v tvare lineárnej kombinácie aproximačných funkcií. V MKP musí približné riešenie na elemente vo všeobecnosti spĺňať nasledujúce požiadavky:
Približné riešenie na elemente je polynóm.
Polynóm musí byť úplný, čo znamená, že musí obsahovať všetky mocniny \(x\) od nultej po najvyššiu.
Derivácie, ktoré sa vyskytujú v slabej formulácii, musia byť vo všeobecnosti pre polynóm nenulové (v prípade slabej formulácie \(\eqref{eq5}\) musí byť nenulová 1. derivácia).
Musia byť splnené interpolačné podmienky pre primárne neznáme.
3.1.1 Lineárne elementy
Najjednoduchšie funkcie, ktoré spĺňajú podmienky 1 až 4 pre slabú formuláciu \(\eqref{eq5}\), sú lineárne polynómy. Približné riešenie na elemente \(\Omega^e\) teda stačí hľadať v tvare
Tieto elementy majú dva uzlové body \(x_1^e = x_A^e\) a \(x_2^e = x_B^e\). Neznáme koeficienty \(c_1, c_2\) určíme tak, aby boli v uzlových bodoch splnené interpolačné podmienky:
s <-3x_nodes <-seq(0, 1, length.out = s +1)# Funkcia na výpočet Lagrangeovho polynómu psi_i(x)lagrange_psi <-function(i, x, x_nodes) { num <-1 den <-1for (j inseq_along(x_nodes)) {if (j != i) { num <- num * (x - x_nodes[j]) den <- den * (x_nodes[i] - x_nodes[j]) } }return(num / den)}# Vykreslenie psi_i(x) pre i = 1, ..., s + 1library(ggplot2)x_plot <-seq(-0.2, 1.2, length.out =200)psi_data <-data.frame(x =rep(x_plot, s +1),y =unlist(lapply(1:(s +1), function(i) sapply(x_plot, lagrange_psi, i = i, x_nodes = x_nodes))),group =factor(rep(1:(s +1), each =length(x_plot))))ggplot(psi_data, aes(x = x, y = y, color = group)) +geom_line(size =1) +geom_vline(xintercept = x_nodes, linetype ="dashed") +geom_hline(yintercept =0, color ="black") +labs(title =paste("Aproximačné funkcie pre s =", s),x ="x", y ="ψ_i(x)", color ="ψ_i") +theme_minimal()
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
3.2 Dosadenie približného riešenia a aproximačných funkcií do slabej formulácie
3.2.1 Lineárny element
Na získanie sústavy rovníc na elemente \(\Omega^e\) použijeme nasledujúci postup.
V slabej formulácii \(\eqref{eq5}\) za neznámu funkciu \(u(x)\) dosadíme približné riešenie \(U^e(x)\) na elemente v tvare \(\eqref{eq8}\)
a za váhovú funkciu \(w(x)\) postupne dosádzame aproximačné funkcie \(\psi_1^e(x)\) a \(\psi_2^e(x)\).
5 Využitie okrajových podmienok pre lineárne elementy
Vo všeobecnosti máme pre úlohu \(\eqref{strong_form}\) zadané okrajové podmienky.
Najprv okrajové podmienky prepíšeme do reči našich neznámych: - v ľavej: \[
U_1 = u(0), \quad Q_1^1 = \left( -a(x) \frac{du}{dx} \right) \bigg|_{x=0}
\] - v pravej: \[
U_{N+1} = u(L), \quad Q_2^N = \left( a(x) \frac{du}{dx} \right) \bigg|_{x=L}
\]
6 Riešenie globálneho systému rovníc a konštrukcia približného riešenia
Na riešenie globálneho systému môžeme použiť: - Priame metódy (napr. Gaussova eliminácia, Thomasov algoritmus), - Iteračné metódy (napr. Gauss-Seidel, SOR, BiCGStab).
Pre lineárne elementy: - Po vyriešení globálneho systému \(\eqref{eq26}\) získame \(U_1, \dots, U_{N+1}\). - Potom na elementoch \(\Omega^e\) s lokálnymi hodnotami \(u_1^e, u_2^e\) priradíme globálne neznáme.
The following object is masked from 'package:ggplot2':
last_plot
The following object is masked from 'package:stats':
filter
The following object is masked from 'package:graphics':
layout
# Parametre úlohyxa <-0# začiatok intervaluxb <-1# koniec intervalun <-4# počet elementovd <-1.0# parameter metódy (implicitný, explicitný, Crank-Nicholson)h <- (xb - xa) / nk <-1.0# parameter na výpočet časového krokuT <-3.0# celkový čas simuláciedeltat <- k * h^2# časový krokNt <-ceiling(T / deltat) # počet časových krokovxVals <-seq(xa, xb, length.out = n +1) # uzly# Funkcia na skladanie globálnych matícvytvor_matrice <-function(n, xVals, a, f) { K <-Matrix(0, nrow = n +1, ncol = n +1) M <-Matrix(0, nrow = n +1, ncol = n +1) F <-numeric(n +1)for (e in1:n) { he <- xVals[e +1] - xVals[e] Ke <-matrix(c(a(xVals[e]) / he, -a(xVals[e]) / he,-a(xVals[e]) / he, a(xVals[e]) / he), 2, 2) Me <-matrix(c(2, 1, 1, 2), 2, 2) * he /6 fe <-f((xVals[e] + xVals[e +1]) /2) * he /2*c(1, 1) idx <-c(e, e +1) K[idx, idx] <- K[idx, idx] + Ke M[idx, idx] <- M[idx, idx] + Me F[idx] <- F[idx] + fe }return(list(K = K, M = M, F = F))}# Definícia funkciía <-function(x) 1+ x # materiálový koeficientf <-function(x) -18* x^2# zdrojový členu0 <-function(x) 1+2* x - x^2# počiatočné podmienky# Okrajové podmienkyd1 <-1.0; B1 <-0.0; g1 <-1.0# ľavý okraj: Dirichletova (ak B1 == 0)d2 <-0.0; B2 <-1.0; g2 <-0.0# pravý okraj: Neumannova (ak d2 == 0)# Inicializácia počiatočného vektoraU <-u0(xVals)# Získanie globálnych matícmatrice <-vytvor_matrice(n, xVals, a, f)K <- matrice$KM <- matrice$MF <- matrice$F# Matrica H pre časový krokH <- M + d * deltat * Kif (B1 ==0) { H[1, ] <-0; H[1, 1] <- d1}if (B2 ==0) { H[n +1, ] <-0; H[n +1, n +1] <- d2}# Ukladanie riešenísolutions <-matrix(0, nrow = Nt, ncol = n +1)solutions[1, ] <- U# Časová slučkafor (t in2:Nt) { G <- (M - (1- d) * deltat * K) %*% U + deltat * Fif (B1 ==0) G[1] <- g1if (B2 ==0) G[n +1] <- g2 U <-as.vector(solve(H, G)) solutions[t, ] <- U}# Príprava údajov na animáciutime_steps <-rep(1:Nt, each = n +1)x_plot <-rep(xVals, Nt)u_plot <-as.vector(t(solutions))plot_data <-data.frame(time = time_steps,x = x_plot,u = u_plot)# Vytvorenie grafufig <-plot_ly() %>%add_lines(x = xVals, y = solutions[1, ], name ="Počiatočný stav", line =list(color ='red', width =3)) %>%add_lines(x = xVals, y = solutions[Nt, ], name ="Konečný stav", line =list(color ='blue', width =3)) %>%add_trace(data = plot_data, x =~x, y =~u, frame =~time,type ='scatter', mode ='lines', name ="Časové kroky",line =list(color ='gray', width =1), showlegend =FALSE) %>%layout(title ="Časový vývoj riešenia",xaxis =list(title ="x"),yaxis =list(title ="u(x, t)"))fig
7 Záver
V tejto práci sme podrobne prešli formuláciou a riešením jednorozmernej okrajovej úlohy pre diferenciálnu rovnicu druhého rádu pomocou metódy konečných prvkov (MKP).Numerické riešenie ukázalo očakávané správanie: pre lineárne aj kvadratické elementy sa približné riešenie rýchlo približuje k presnému riešeniu. Animácia časového priebehu pomohla vizualizovať dynamiku problému a overiť stabilitu metódy.
Práca poskytuje pevný základ na rozšírenie riešenia na viacrozmerné problémy alebo úlohy vyššieho rádu, prípadne na využitie pokročilejších aproximačných funkcií a nelineárnych metód. (1, 2)
