FEM in SciLab

Aim: Temperature distribution FEM simulation in SciLab.

Temperature distribution has been defined with different heat sources and dissipators. Point matrix has been implemented into SciLab. Project has been done step-by-step, what has resulted in working program.

//Analizowany obszar

SA=12; //Szerokosc obszaru [m] SB=8; //Wysokosc obszaru [m] D1=0.1; //Gestosc siatki – dzielnik SA i SB
D2=0.2; //Odleglosc punktow zrodlowych od obszaru

N1=189;
N2=211;
N3=0;

NW=(SA/D1-1)*(SB/D1-1); //Zapis macierzy o rozmiarze NW ( liczba puktow, w ktorych obliczamy temperature )

xw=zeros(NW); //Wspolrzedne x punktow, w ktorych obliczamy temperature
k=1;
for i=1:(SB/D1-1)
for j=1:(SA/D1-1)
xw(k)=1+D2+j*D1;
k=k+1;
end
end

yw=zeros(NW); //Wspolrzedne y punktow, w ktorych obliczamy temperature
k=1;
for i=1:(SA/D1-1)
for j=1:(SB/D1-1)
yw(k)=1+D2+j*D1;
k=k+1;
end
end

T=zeros(NW); //Macierz wynikowa

N=2*(SA/D1+SB/D1); //Liczba rownan

//Zapis macierzy o rozmiarze N ( liczba punktow zrodlowych )

xz=zeros(N); //Wspolrzedne x punktow zrodlowych
xz(1)=1;
for i=2:(SA/D1)
xz(i)=xw(i-1);
end
for i=(SA/D1+1):(SA/D1+1+SB/D1)
xz(i)=1+SA+2*D2;
end
for i=(SA/D1+2+SB/D1):(2*SA/D1+1+SB/D1)
xz(i)=xz(2*SA/D1-i+2+SB/D1);
end
for i=(2*SA/D1+2+SB/D1):(2*SA/D1+2*SB/D1)
xz(i)=1;
end

yz=zeros(N); //Wspolrzedne y punktow zrodlowych
for i=1:(SA/D1+1)
yz(i)=1;
end
for i=(SA/D1+2):(SA/D1+SB/D1)
yz(i)=yw(i-(SA/D1+1));
end
for i=(SA/D1+SB/D1+1):(2*SA/D1+SB/D1+1)
yz(i)=yz(SA/D1+SB/D1)+D1+D2;
end
for i=(2*SA/D1+SB/D1+2):(2*SA/D1+2*SB/D1)
yz(i)=yz(SA/D1+SB/D1-i+2*SA/D1+SB/D1+2);
end

//Zapis macierzy o rozmiarze N1 ( liczba punktow kollokacji na brzegu g1 – war. Dirichleta – znana temperatura na brzegu )

xc1=zeros(N1); //Wspolrzedne x punktow kolokacji
xc1(1)=1+D2;
for i=2:(SA/D1-11)
xc1(i)=xc1(1)+(i+10)*D1;
end
for i=(SA/D1-10):(SA/D1-11+SB/D1)
xc1(i)=1+D2;
end

yc1=zeros(N1); //Wspolrzedne y punktow kolokacji
yc1(1)=1+D2;
for i=2:(SA/D1-11)
yc1(i)=1+D2;
end
yc1(SA/D1-10)=1+D2+SB;
for i=(SA/D1-9):(SA/D1-11+SB/D1)
yc1(i)=yc1(i-1)-D1;
end

Tb=zeros(N1); //Wartosc temperatury w punktach
Tb(1)=15;
for i=2:(SA/D1-11)
Tb(i)=Tb(1)+(i-2)*(5/(SA/D1-13));
end
for i=(SA/D1-10):(SA/D1-11+SB/D1)
Tb(i)=20-(i-(SA/D1-10))*(5/((SA/D1-11+SB/D1)-(SA/D1-10)));
end

//Zapis macierzy o rozmiarze N2 ( liczba punktow kollokacji na brzegu g2 – war. Neumana – znany jest strumien ciepla )

xc2=zeros(N2); //Wspolrzedne x punktow kolokacji
for i=1:11
xc2(i)=1+D2+i*D1;
end
for i=12:92
xc2(i)=1+D2+SA;
end
for i=93:211
xc2(i)=xc2(i-1)-D1;
end

yc2=zeros(N2); //Wspolrzedne y punktow kolokacji
for i=1:12
yc2(i)=1+D2;
end
for i=13:92
yc2(i)=yc2(i-1)+D1;
end
for i=93:211
yc2(i)=yc2(92);
end

nx2=zeros(N2); //Skladowe x jednostkowego wektora prostopadlego do brzegu)
for i=1:11
nx2(i)=xc2(i);
end
for i=12:92
nx2(i)=1+SA+1.5*D2;
end
for i=93:211
nx2(i)=xc2(i);
end

ny2=zeros(N2); //Skladowe y jednostkowego wektora prostopadlego do brzegu)
for i=1:11
ny2(i)=1+0.5*D2;
end
for i=12:92
ny2(i)=yc2(i);
end
for i=93:211
ny2(i)=1+1.5*D2+SB;
end

qn=zeros(N2); //Wartosc strumienia ciepla w punktach
for i=1:11
qn(i)=10;
end
for i=12:92
qn(i)=0;
end
for i=93:211
qn(i)=-15;
end
for j=1:6
for i=1:10
qn(102+20*(j-1)+i)=4;
end
end

//Zapis macierzy o rozmiarze N3 ( liczba punktow kollokacji na brzegu g3 – war. Robina – strumien ciepla proporcjonalny do roznicy temperatury na brzegu i otoczenia )

xc3=zeros(N3); //Wspolrzedne x punktow kolokacji
yc3=zeros(N3); //Wspolrzedne y punktow kolokacji
nx3=zeros(N3); //Skladowe x jednostkowego wektora prostopadlego do brzegu)
ny3=zeros(N3); //Skladowe y jednostkowego wektora prostopadlego do brzegu)
Bi=zeros(N3); //Wartosc liczby Biota

A=zeros(N,N); //Macierz wspolczynnikow
B=zeros(N); //Wektor prawej strony ukladu rownan
X=zeros(N); //Wektor niewiadomych

for i=1:N1 //Kolokacja warunku na brzegu g1
B(i)=Tb(i);
for j=1:N
A(i,j)=log((xc1(i)-xz(j))^2+(yc1(i)-yz(j))^2);
end
end

for i=1:N2 //Kolokacja warunku na brzegu g2
B(N1+i)=qn(i);
for j=1:N
A(N1+i,j)=2*(xc2(i)-xz(j))/((xc2(i)-xz(j))^2+(yc2(i)-yz(j))^2)*nx2(i)+2*(yc2(i)-yz(j))/((yc2(i)-yz(j))^2+(xc2(i)-xz(j))^2)*ny2(i);
end
end

for i=1:N3 //kolokacja warunku na brzegu g3
B(N1+N2+i)=0.0;
for j=1:N
A(N1+N2+i,j)=2*(xc3(i)-xz(j))/((xc3(i)-xz(j))^2+(yc3(i)-yz(j))^2)*nx3(i)+2*(yc3(i)-yz(j))/((yc3(i)-yz(j))^2+(xc3(i)-xz(j))^2)*ny3(i)+Bi(i)*log((xc3(i)-xz(j))^2+(yc3(i)-yz(j))^2);
end
end

X=linsolve(A,-B); //Rozwiazywanie ukladu rownan liniowych

Y=zeros(NW); //Odkrecenie Y na normalna postac
k=1;
for i=1:(SB/D1-1)
for j=1:(SA/D1-1)
Y(k)=1+D2+SB-i*D1;
k=k+1;
end
end

for k=1:NW //Tworzenie macierzy rozwiazan
W=0;
for i=1:N
W=W+X(i)*log((xw(k)-xz(i))^2+(Y(k)-yz(i))^2);
end
T(k)=W;
end

k=1;
SRF=zeros((SA/D1-1),(SB/D1-1)); //Macierz temperatur
for i=1:(SB/D1-1)
for j=1:(SA/D1-1)
SRF(j,(SB/D1-i))=T(k);
k=k+1;
end
end

n1=zeros(SA/D1-1); //Macierz wspolrzednych y
for i=1:(SA/D1-1)
n1(i)=1+D2+i*D1;
end

n2=zeros(SB/D1-1); //Macierz wspolrzednych x
for i=1:(SB/D1-1)
n2(i)=1+D2+i*D1;
end

f=scf(0); //Plot wynikow
plot3d1(n1,n2,SRF)
f.color_map = jetcolormap(1000);