Transcript Sveučilište u Zagrebu Fakultet kemijskog inženjerstva i tehnologije Kolegij: Matematičke metode u kemijskom inženjerstvu METODE KONAČNIH ELENMENATA I KONAČNIH RAZLIKA Neven Ukrainczyk [email protected] Zagreb, travanj 2003. I BAZNE.

Sveučilište u Zagrebu
Fakultet kemijskog inženjerstva i tehnologije
Kolegij: Matematičke metode u kemijskom inženjerstvu
METODE KONAČNIH ELENMENATA I
KONAČNIH RAZLIKA
Neven Ukrainczyk
[email protected]
Zagreb, travanj 2003.
I BAZNE FUNKCIJE KONAČNIH ELEMENATA
a)
b)
Distribucija temperature u(x) po duljini šipke.
Polinomska funkcijska veza nađena metodom najmanjih kvadrata
– neprihvatljive oscilacije između točaka.
u(x) = a + b x + c x2+ d x3+...,
a) Izmjerena temperatura u u ovisnosti o duljini luka s.
b) Podjela domene na tri elementa u kojima linearni
polinomi opisuju ovisnost.
Linearne bazne funkcije
-zamjena parametara s vrijednostima funkcije
u na granicama elemenata:
u( )  (1   ) u1   u 2
gdje je (0  1) normalizirana mjera udaljenosti na krivulji.
1 ( ) 1  
 2 ( )  
u( ) 1 ( ) u1   2 ( ) u2
Linearne bazne funkcije
1 ( ) 1  
 2 ( )  
Odnos globalnih i lokalnih čvorova.
u n  U  ( n ,e )
u( ) 1 ( ) u1   2 ( ) u2
Izmjereno temperaturno polje opisano čvornim parametrima i
linearnim baznim funkcijama, koje je sad neprekinuto na spojištima
elemenata.
Bazne funkcije – težinske funkcije
Bazne funkcije se mogu smatrati kao težinske
funkcije čvornog parametra:
0
1

4
1

2
3

4
 1
u(0)  (1  0) u1  0 u2
1
1
1
3
1
u ( )  (1  ) u1  u 2  u1  u 2
4
4
4
4
4
1
1
1
1
1
u ( )  (1  ) u1  u 2  u1  u 2
2
2
2
2
2
3
3
3
1
3
u ( )  (1  ) u1  u 2  u1  u 2
4
4
4
4
4
u(1)  (1  1) u1  1u 2  u 2
(a)...(b) Težinske funkcije n pridodijeljene
globalnim čvorovima n = 1...4.
u( ) 1 ( ) u1   2 ( ) u2
x( ) 1 ( ) x1   2 ( ) x2
u ( )    n ( ) u n
n
x( )    n ( ) x n
n
Veza u i x preko normalizirane kordinate elementa .
Kvadratna bazna funkcija
za kvadratnu ovisnost u na elementu su potrebna tri
čvorna parametra u1, u2 i u3 :
u( ) 1 ( ) u1   2 ( ) u2  3 ( ) u3
Jednodimenzionalne kvadratne bazne funkcije
Dvodimenzionalna bilinearna bazna funkcija (linearna na
1 i 2 koordinati) je konstruirana od produkta prije
navedenih jednodimenzionalnih linearnih funkcija:
u(1 , 2 ) 1 (1 , 2 ) u1  2 (1 , 2 ) u2  3 (1 ,  2 ) u3  4 (1 , 2 ) u4
gdje je
1 (1 ,  2 )  (1  1 ) (1   2 )
 2 (1 ,  2 )  1 (1   2 )
 3 (1 ,  2 )  (1  1 )  2
 4 (1 ,  2 )  1  2
Dvodimenzionalne bilinearne bazne funkcije.
II STACIONARNO VOĐENJE TOPLINE
Jednodimenzionalno stacionarno vođenje topline
Iz jednostavne bilance topline infinitezimalnog dijela materijala
dobivamo:
Promjena toplinskog fluksa = generacija topline
d
(toplinski fluks)  gubitak topline  0
dx
d  du 
 k
  q(u, x)  0
dx  dx 
gdje je u temperatura, x duljina štapa, gubitak topline i k toplinska
vodljivost (W/(moC)).
Razmotrimo slučaj q = u
d  du 
 k
u  0
dx  dx 
0  x 1
s graničnim uvjetima: u(0) = 0 i u(1) = 1.
Ova jednadžba (uz k = 1) ima egzaktno rješenje
e
u ( x)  2
(e x  e  x )
e 1
1.
2.
3.
4.
5.
6.
7.
8.
Da bi riješili jednadžbu metodom konačnih
elemenata potrebni su ovi koraci:
Pisanje jednadžbe u integralnom obliku.
Parcijalna integracija (1D) ili korištenje Greenovog
teorema (2D i 3D) za snižavanje reda derivacije.
Aproksimacija temperaturnog polja konačnim
elementima.
Integracije na elementima za izračunavanje elementna
matrice toplinske vodljivosti i vektora toplinskog
toka.
Slaganje globalne jednadžbe.
Primjena graničnih uvjeta.
Rješenje globalne jednadžbe.
Izračunavanje toplinskih tokova.
1. Integralna jednadžba
 R  dx  0
gdje je R ostatak
d  du 
R  k   u
dx  dx 
 d  du 

0  dx  k dx   u  dx  0
1
2. Parcijalna intergracija
Supstitucijom u= i
parcijalnu integraciju:
u
du
v  k
dx
u formulu za
dv
du
dx  u v   v
dx
dx
dx
1
 
d 
du 
du  
du d 



k
dx



k


k
0 dx  dx    dx   0  dx dx dx
0
1
1
1
  du d 

 du 

k

u


dx


0   dx dx 
k dx  

0
1
3. Aproksimacija konačnim elementima
Podjelimo domenu 0 < x < 1 na tri jednaka elementa
i zamjenimo kontinuiranu veličinu u(x) na svakom
pojedinom elementu parametarski zadanom
aproksimacijom konačnih elemenata:
u( )  1 ( ) u1  2 ( ) u2  n ( ) un
n
x( ) 1 ( ) x1  2 ( ) x2  n ( ) xn
n
gdje su
1 ( ) 1   su  2 ( )  
Za test funkciju biramo  = m (Galjerkinova
pretpostavka)
To prisiljava da ostatak R bude okomit na bazne
funkcije
1
2
1
3
3
1
0
0
1
3
2
3
 ..dx   ..dx   ..dx   ..dx
x2
1
x1
0
..
dx

..
J
d



J
dx
d
gdje je
Jacobijeva determinanta za
transformaciju iz x-koordinata u -koordinate.
4. Integracije na elementima
 du d

0  k dx dx  u   J d
1
gdje je
u  n un i  m .
n
 d m d  d n d 

n un   k d dx d dx  nm  J d

0
1
Iz odnosa, x: = 1:3, Jacobian je
dx 1

d 3
1
 d m d d n d

 d m d n

   k
  n m  J d    k
3
3   n m  J d
d dx d dx
d
d


0
0
1
Emn
J
Da izračunamo Emn, uvrstimo bazne funkcije:
d1
ili
 1
d
d 2
ili
1
d
1 ( )  1  
 2 ( )  
2
1


1   d1 
1
1
1
2
2
2
E11   9 k    1 d   9 k (1)  (1   ) d   9 k  

3 0   d 
3
3
3


0

1


1
1
E12  E21   9 k  
3
6
1
1
E22   9 k  
3
3
Emn
1 
1
3 9 k  3 




1 
1
  9 k  
6
3 
1
1 
 9 k   
3
6 
1
1 
9k   
3
3 
5. Slaganje globalne jednadžbe
Redovi globalne matrice su formirani od globalnih
težišnih funkcija
Slaganjem dobivamo (k = 1):
 28
 9

  53
 18

 0


 0

53

18
28 28

9
9
53

18
0
53

18
28 28

9
9
53

18
0

0 
 U 
1

0  
 U 2 
   = Vektor toplinskog
53  U 3
toka

 
18  U 4 
28 

9 
1
  du d 

 du 
0   k dx dx   u  dx  k dx  
0
1
x 1
 du 
 du 
 du 
k


k





k
 dx 

 x 0  dx  x 1  dx  x 0
1 x0  1
1 x1  0
x 1
 du 
 du 
 ulaz toplinskog toka za čvor 1
k dx 1    k dx 

 x 0

 x 0
Slično, za čvor 2 i 3:
x 1
x 1
 du 
k dx  n   0
x 0
(čvor 2 i 3)
 du 
 du 
 k dx 4    k dx   ulaz toplinskog toka za čvor 4.

 x 0 
 x 1
Sastavljanjem tih globalnih jednadžbi dobiva se
 28
 9

  53
 18

 0


 0

ili
53

18
28 28

9
9
53

18
0
0
53

18
28 28

9
9
53

18

0 
  du  

k


 
 U 
 dx  x 0 
1

0   

0
 U 2 

   
53  U 3
0



 

18  U 4    du 
 k



28
  dx  x 1 

9 
Ku f
6. Granični uvjeti
Granični uvjeti su primijenjeni direktno na prvi i
zadnji čvor
U1
0
53
56
53
 U1  U 2  U 3
0
18
9
18
53
56
53
 U 2  U3  U 4  0
18
9
18
U 4 1
7. Rješenje
Rješenje sustava jednadžbi je:
U2 = 0.2885
U2 egz= 0.2882
U3 = 0.6098
U2 egz = (0.6102)
8. Toplinski tokovi
Toplinski tokovi za čvor 1 i 4 su izračunati uvrštenjem
čvornih rješenja U1 = 0, U2 = 0.2885, U3 = 0.6098 i
U4 = 1
Rješenje 1D vođenja topline metodom konačnih elemenata.
III NEUSTALJENO VOĐENJE TOPLINE
Metoda konačnih razlika
Eksplicitni oblik
Neustaljeno jednodimenzionalno vođenje topline:
u
 2u
k 2
t
x
(0  x  L, t  0)
gdje je k toplinska vodljivost, a u = u(x,t) temperatura,
uz granične uvjete U(0,t) = u0 i u(L,t) = u1
i početne uvjete u(x,0) = 0
Mreža konačnih razlika za rješavanje nestacionarnog 1D vođenja topline.
Jednadžba je centrirana na čvor mreže (i,n) označen kružićem. Područje
već poznatog rješenja do n-tog koraka je lagano osjenćano.
u( x, t )  u(i x, n t )  u
n
i
u
n
i 1
2
3
n
3
n
n
i
n
i 1
n
1 2  u 
1 3  u 
 u 
 u  x   x  2   x  3   O(x 4 )
 x  i 2
 x  i 6
 x  i
n
u
n
1 2  u 
1 3  u 
 u 
 u  x   x  2   x  3   O(x 4 )
 x  i 2
 x  i 6
 x  i
n
2
n
i
 u 
 u  t    O(t 2 )
 t  i
n
u
n 1
i
n
i
gdje O(x ) i O(t ) predstavlja sve ostale izraze
u Taylorovom redu.
4
2
u
n
i 1
u
n
i 1
n
 u
 2 u  x  2   O(x 4 )
 x  i
2
n
i
2
n
n
n
n
 u
ui 1  2 ui  ui 1
2


O
(

x
)
 2
2
x
 x i
2
n 1
i
 u  u  u
 O(t )
  
t
 t i
n
n
i
n
n
n
u in 1  u in
u

2
u

u
4
i
i 1
 O(t 2 )  k i 1

O
(

x
)
2
t
x
u
n 1
i
t
 u  k 2 (u in1  2 u in  u in1 )  O(t 2 , x 2 )
x
n
i
Stabilnost rješenja
t
1
k 2 
x
2
2
(1) n2 2t
u( x, t )  x  
e
sin(nx)
 n1 n

n
Analitička i numerička rješenja nestacionarnog 1D
vođenja topline prikazuju efekt veličine x i t.
Nestacionarna advekcijsko-difuzijska jednadžba
prijenosa:
u
 v u  k 2 u  f
t
gdje je u temperatura (ili koncentracija), advektivni
transport (u slučaju konvekcijskog prijenosa topline)
sa poljem brzina v, k je koeficijent toplinske
vodljivosti (ili koeficijent difuzije) a f je izraz za
izvor (generaciju).
u( x, t ) 
M
4t
e
( x vt ) 2
4 Dt
Advekcio-difuzijski odziv na jedinični impuls. x = 0.1, t
= 0.001 s za 0<t<0.01 s i t = 0.01 s za t  0.01 s.