cieplo.m

function cieplo()
 
 
%% SKRYPT OCTAVE / CALFEM DO ANALIZY MES USTALONEGO PRZEPŁYWU CIEPŁA
% AUTOR: Piotr Pluciński, PK WIL L-10
 
clf
 
Edof=[ 1 1 2 3 ; 2 1 3 4 ] ; % nr_elementu, numery_stopni_swobody_w_elementach
 
K=zeros(4); 
P=zeros(4,1);
kx = 4; % współczynnik przewodności cieplnej wzdłuż x [J / Cms]
ky = 4; % współczynnik przewodności cieplnej wzdłuż y [J / Cms]
f = 45; % intensywność generacji ciepła wewnątrz ciała [J / m2s]
D=[kx 0 ; 0 ky]; 
 
Coord=[0 0 ; 2 0 ; 2 1 ; 0 1]; % współrzędne węzłów
 
Dof=(1:4)';  % stopnie swobody globalne w poszczególnych węzłach(1:ilość węzłów)
 
[Ex,Ey]=coordxtr(Edof,Coord,Dof,3); 
 
 
figure(1)
eldraw2(Ex,Ey, [1,2,2], Edof(:,1)); 
hold on;
grid on;
axis equal;
title('DYSKRETYZACJA MES ZADANIA USTALONEGO PRZEPŁYWU CIEPŁA');
xlabel('x');
ylabel('y');
 
%% 
 
%[K1,P1]=flw2te(Ex(1,:),Ey(1,:),1,D,f);
%[K2,P2]=flw2te(Ex(2,:),Ey(2,:),1,D,f);
 
%[K,P]=assem(Edof(1,:),K,K1,P,P1);
%[K,P]=assem(Edof(2,:),K,K2,P,P2);
 
 for el=1:size(Edof,1)
     [Ke,Pe]=flw2te(Ex(el,:),Ey(el,:),1,D,f);
     [K,P]=assem(Edof(el,:),K,Ke,P,Pe);
 end
 
% budowa wektora prawej strony układu równań MES
P=P+flux(Coord,[3 4],-30); %% 4 nr stopni swobody na brzegu ze definiowanym strumieniem, -30 wartość strumienia,
% znakowanie strumienia: dodatni gdy ciepło wpływa do ciała; ujemny gdy ciepło wypływa z ciała
%% ile mamy brzegów ze strumieniem tyle razy trzeba powtórzyć powyższą linię
 
bc=[ 2 10; 3 10]; % wektor warunków brzegowych (nr węzła/ss - wartość temp.)
 
%% ROZWIĄZANIE UKŁADU RÓWNAŃ MES i RYSUNEK TEMPERATURY
[T,R]=solveq(K,P,bc) 
 
 
 
%% POWRÓT DO ELEMENTU - WYZNACZENIE WARTOŚCI WEKTORA STRUMIENIA CIEPŁA
Te=extract_ed(Edof,T);
 
%Q(1,:)=flw2ts(Ex(1,:),Ey(1,:),D,Te(1,:));
%Q(2,:)=flw2ts(Ex(2,:),Ey(2,:),D,Te(2,:));
 
 
 for el=1:size(Edof,1)
     Q(el,:)=flw2ts(Ex(el,:),Ey(el,:),D,Te(el,:));
 end
%% 
 
Q
 
% mapa temperatur 
sfac=scalfact2(Ex,Ey,Q,0.25);
 
figure(2)
filloct(Ex',Ey',Te')
title(['TEMPERATURE [^oC] : T_m_a_x = ' num2str(max(T)) ', T_m_i_n = ' num2str(min(T))]);
xlabel('x');
ylabel('y');
zlabel('T');
colorbar
colormap jet
axis equal
 
 
figure(3)
eldraw2(Ex,Ey, [1,2,2]); 
elflux2(Ex,Ey,Q,[1,4],sfac);
eliso2(Ex,Ey,Te,10,[2,1]);
axis equal
 
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
function Z=flux(coord,wez,q)
 
Z=zeros(size(coord,1),1)
L=sqrt((coord(wez(1),1)-coord(wez(2),1))^2+(coord(wez(1),2)-coord(wez(2),2))^2)
Z(wez)=q*L/2;