www.pudn.com > 87361026FEA.rar > fem_example2.m, change:2006-02-28,size:1591b


% fem_example.m
%
% Jack Chessa 
% Northwestern University
%
% A simple 1D finite element program in Matlab
%
%  
%    hint( Tint - T1 ) =============================     hext( Text - Tnn )
%        Tint  <~~~~~~ |          k A              | ~~~~~~~>      Text                           
%                      =============================
%  

%--------------------------- INPUT SECTION -------------------------------%
% define mesh
node=[0.0 1.5 7.2 7.5 8.0 9.0 10.0]'; % node coordinates
element=[1 2;2 3;3 4;4 5;5 6;6 7];    % element connectivity

% define problem parameters
k=100*[1 8 1 2 2 2];  % element thermal conductivities
a=    [1 1 1 1 1 1];  % element cross sectional areas

Text=0;     % temperature outside
Tint=100;   % temperature inside

hext = 1;   % conductivity coefficeint outside
hint = 10;   % conductivity coefficeint inside

%------------------------- END OF INPUT SECTION --------------------------%
nn=size(node,1);    % number of nodes
ne=size(element,1); % number of elements

ka=k.*a;            % element area conductivities

% define fe operators
K=zeros(nn,nn);        % thermal stiffness matrix
f=zeros(nn,1);         % force vector
d=zeros(nn,1);         % nodal displacement

% compute stiffness matrix
for e=1:ne
  
  sctr = element(e,:);  
  
  le=abs( node(sctr(2)) - node(sctr(1)) ); 
  ke=ka(e)/le*[1 -1;-1 1];
  
  K(sctr,sctr) = K(sctr,sctr) + ke;
  
end

% thermal load ( convection conditions )
f(nn)= -hext*Text;
K(nn,nn)=K(nn,nn)-hext;
f(1)= -hint*Tint;
K(1,1)=K(1,1)-hint;

% solve system
d=K\f;

% plot results
clf
plot(node,d,'r-o')