Przykład 1 modelowania jednowymiarowego przepływu ciepła 1.

1. Model przepływu przez ścianę wielowarstwową

Ściana składa się trzech warstw o różnych grubościach wykonana z różnych materiałów. Na jednej ze ścian zewnętrznych dana jest temperatura T0 na przeciwległej znana jest temperatura otoczenia Tot

oraz współczynnik wymiany ciepła α. Określić temperatury w kolejnych warstwach ściany:

Dane:

λ₁ = 20 W/m °C, λ2 = 30 W/m °C, λ3 = 50 W/m °C, L1=0.2m, L2=0.1m, L1=0.1m α = 25 W/m² °C, Tot = 800 °C, T0 = 20 °C

2. Dyskretyzacja - podział na elementy skończone Wyróżniono 3 liniowe elementy skończone (4 węzły)

3. Wyznaczenie lokalnych macierzy przewodności cieplnej

Lokalna macierz przewodności cieplnej dla liniowego elementu 1D:

gdzie: λ_e – przewodność cieplna elementu skończonego, , Le – długość elementu skończonego

W przykładzie należy wyznaczyć 3 lokalne macierze przewodności cieplnej , ponieważ wyróżniono 3 elementy skończone

1 1

1 1

e e

Le

λ  − 

=   −  

1 1 2 2 3 3

1 1 11 12 2 2 11 12 3 3 11 12

1 1 2 2 3 3

1 21 22 2 21 22 3 21 22

1 1 1 1 1 1

... ...

1 1 1 1 1 1

k k k k k k

L k k L k k L k k

λ

λ  −    λ  −     −   

=   −   =     =   −   =     =   −   =    

k k k

4. Agregacja

Zagregowana macierz przewodności cieplnej dla przykładu ma postać:

5. Uwzględnienie warunków brzegowych – metoda kary

Zgodnie ze schematem metody kary dla zadania przepływu ciepła (str.2) modyfikacji podlegają elementy leżące na przekątnej globalnej macierzy przewodności cieplnej oraz prawa strona równania (macierz strumienia ciepła R). W węźle nr 1 występuje warunek brzegowy 3 rodzaju (WB3), natomiast w węźle nr 4 warunek brzegowy 1pierwszego rodzaju (WB2):

6. Rozwiązanie układu równań MES

Zastosowanie metody kary umożliwia potraktowanie wektora temperatur węzłowych jako niewiadomej w równaniu MES:

Odpowiada to przekształceniu równania do postaci

1 1

11 12

1 1 2 2

21 22 11 12

2 2 3 3

21 22 11 12

3 3

21 22

0 0

0 0

0 0

k k

k k k k

k k k k

k k

 

 + 

 

=  + 

 

 

 

K

1 1

11 12 1

1 1 2 2

21 22 11 12 2

2 2 3 3

21 22 11 12 3

3 3

0

21 22 4

0 0

0 0 0 0

0 0

T

k k T

k k k k T

k k k k T

C T

k k C T

α α

 +     

 +     

    =  

 +     

     

+    

     

 

KT = R

T = K R-1

7. Kod do programu SciLab dla powyższego przykładu

clear

clc //określenie długości poszczególnych elementów

L1=0.2;

L2=0.1;

L3=0.1;

//określenie współczynników przewodności cieplnej lam1=20;

lam2=30;

lam3=50;

//określenie warunków brzegowych alfa=25;

Tot=800;

T0=20;

//wyznaczenie macierzy przewodności cieplnej poszczególnych elementów K1=lam1/L1*[1, -1; -1, 1]

K2=lam2/L2*[1, -1; -1, 1]

K3=lam3/L3*[1, -1; -1, 1]

//Zainicjowanie globalnej macierzy sztywności KG=zeros(4,4)

//Agregacja globalnej macierzy sztywności KG(1:2,1:2)=KG(1:2,1:2)+K1

KG(2:3,2:3)=KG(2:3,2:3)+K2 KG(3:4,3:4)=KG(3:4,3:4)+K3

//Zainicjowanie macierzy kolumnowej węzłowych strumieni ciepła R=zeros(4,1)

//Uwzględnienie warunków brzegowych - metoda kary C=max(KG)*10^4; //Określenie stałej C

R(1,1)=R(1,1)+(alfa*Tot) R(4,1)=R(4,1)+(C*T0) KG(1,1)=KG(1,1)+alfa KG(4,4)=KG(4,4)+C //Rozwiązanie układu równań T=inv(KG)*R;

//wyświetlenie wektora węzłowych temperatur T disp(T)

Przykład 2 modelowania jednowymiarowego przepływu ciepła (uwzględnienie wewnętrznych źródeł ciepła)

1. Model jednowymiarowego przepływu ciepła z wewnętrznym źródłem ciepła.

W płycie o grubości 25cm i współczynniku przewodności cieplnej λ=0.8W/m °C generowane jest ciepło w wydajności 4000 W/m³. Zewnętrzną powierzchnię płyty z obydwu stron opływa powietrze o temperaturze 30°C. Współczynnik wnikania ciepła wynosi α=20W//m² °C. Wyznaczyć rozkład temperatury wewnątrz płyty.

Dane:

λ = 0.8 W/m °C, α = 20 W/m² °C, Tot = 30 °C, Q = 4000 W/m³ 2. Dyskretyzacja - podział na elementy skończone

Zadanie to można zamodelować jako zadanie symetryczne, ponieważ oś pionowa prostopadła do płyty przechodząca przez jej środek dzieli układ na dwie symetryczne połówki. Symetria dotyczy geometrii oraz warunków brzegowych. W zagadnieniach przepływu ciepła w węzłach leżących na osi symetrii zadaje się warunek adiabatyczny (izolacji), tj. warunek brzegowy drugiego rodzaju w którym strumień ciepła q jest równy 0. Symetryczną połowę układu podzielono na 2 elementy skończone.

3. Wyznaczenie lokalnych macierzy przewodności cieplnej

Lokalna macierz przewodności cieplnej dla obydwu elementów będzie identyczna ponieważ elementy posiadają ten sam współczynnik przewodności cieplnej oraz tę samą długość, tj. 0,0625m:

4. Agregacja

Zagregowana macierz przewodności cieplnej dla przykładu ma postać:

5. Wyznaczenie oraz zagregowanie macierzy kolumnowej strumienia ciepła

Ponieważ wydajność wewnętrznego źródła ciepłą jest równa 4000 W/m, to w każdym elemencie o długości 0.0625m wydzieli się 250W ciepła (Q*L=4000*0.0625). Każdy z zastosowanych tu elementów skończonych posiada 2 węzły więc w każdym z węzłów należy zadać 125W ciepła.

Agregacja powoduje, że dla węzła wspólnego elementom skończonym 1 i 2 wartość strumienia ciepła w węźle wynosi 250W.

1,2

1 1

1 1

λ  − 

=   −  

1 1

11 12

1 1 2 2

21 22 11 12

2 2

21 22

0

0

k k

k k k k

k k

 

 

=  + 

 

 

K

6. Uwzględnienie warunków brzegowych – metoda kary Zgodnie ze schematem metody kary dla zadania przepływu ciepła

- w globalnej macierzy przewodności cieplnej K do elementu o indeksie (3,3) dodawany jest współczynnik α

- do globalnego wektora strumienia ciepła R do elementu o indeksie (3,1) dodawany jest iloczyn α i Tot

8. Rozwiązanie układu równań MES

Podobnie jak w poprzednim przykładzie rozwiązujemy układ równań:

z którego wyznaczamy węzłowe wartości temperatur

9. Kod do programu SciLab dla powyższego przykładu

clear

clc //określenie długości poszczególnych elementów

L1=0.0625;

L2=0.0625;

//określenie współczynników przewodności cieplnej lam1=0.8;

lam2=0.8;

//określenie warunków brzegowych alfa=20;

Tot=30;

//zdefiniowanie wydajności wewnętrznych źdódeł ciepła Q=4000;

//wyznaczenie macierzy przewodności cieplnej poszczególnych elementów K1=lam1/L1*[1, -1; -1, 1]

K2=lam2/L2*[1, -1; -1, 1]

//Zainicjowanie globalnej macierzy sztywności KG=zeros(3,3)

//Agregacja globalnej macierzy sztywności KG(1:2,1:2)=KG(1:2,1:2)+K1

KG(2:3,2:3)=KG(2:3,2:3)+K2

//Zainicjowanie macierzy kolumnowej węzłowych strumieni ciepła globalnej oraz lokalnych dla elementów R=zeros(3,1)

R1=zeros(2,1) R2=zeros(2,1)

//Uwzględnienie wewnętrznego źródła ciepła w macierzy kolumnowej węzłowych strumieni ciepła R1(1:2,1)=Q*L1/2

R2(1:2,1)=Q*L2/2

//Agregacja globalnej macierzy strumeni ciepła R(1:2,1)=R(1:2,1)+R1

R(2:3,1)=R(2:3,1)+R2

//Uwzględnienie warunków brzegowych - metoda kary C=max(KG)*10^4; //Określenie stałej C

R(3,1)=R(3,1)+(alfa*Tot) KG(3,3)=KG(3,3)+alfa //Rozwiązanie układu równań T=inv(KG)*R;

//wyświetlenie wektora przemieszczeń Q disp(T)

/ 2 0 / 2 125

/ 2 / 2 250

0 / 2 / 2 125

QL QL

R QL QL QL

QL QL

       

       

=    +   =   = 

       

       

KT = R