%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% One dimensional finite element implementation for steady state heat
% solve using the problem given.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [T,reaction] = SteadyHeatForStudents()
clear % clear the workspace variables
clc % clear the command window
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Problem data
%
% length of ``bar''
%
Length = 4.0;
%
% cross sectional area (possibly a function of position)
%
Area = @(x) 0.1;
%
% thermal conductivity (possibly a function of position)
%
Ktherm = @(x) 2.0;
%
% heat source/sink (possibly a function of position)
%
source = @(x) 5.0;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Create the mesh you want
%
% number of elements
%
nElem = 2;
%
% nodes per element
%
nNode = 2;
%
% number of integration points per element
%
nInt = 1;
%
% the total number of nodes in this mesh
%
nNodes = nElem*nNode - (nElem - 1);
%
% generate the mesh
%
coord = zeros(nNodes,1);
for i=2:nNodes;
coord(i) = coord(i-1) + Length/(nNodes-1);
end
%
% here is the element to nodal connectivity in the form
% node on left = connect(elem,1)
% next node moving right = connect(elem,2)
% ....
% node on the right = connect(elem,nNode)
%
connect = zeros(nElem,nNode);
for i=1:nNode;
connect(1,i) = i;
end
for i=2:nElem;
for j=1:nNode;
connect(i,j) = connect(i-1,j) + nNode - 1;
end
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% init
%
nodeCoords = zeros(nNode,1);
K = zeros(nNodes,nNodes);
F = zeros(nNodes,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% initial condition on the temperature (not needed for steady state)
%
%T = 0.0*ones(nNodes,1);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Boundary conditions
%
% flux boundary condition
%
flux = -5.0; % magnitude AND direction of flux BC
elemFlux = nElem; % element with the flux BC
nodeFlux = nNode; % local node number to apply the flux BC
%
% temperature boundary condition
%
Tbc = 0.0;
nodeT = 1; % global node number to apply the temperature BC
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% loop over elements
%
for elem=1:nElem;
% get the coordinates for the nodes on this element
%
for i=1:nNode;
nodeCoords(i) = coord(connect(elem,i));
end
% compute this elements stiffness and force
%
[Ke,Fe] = element(nInt,nNode,nodeCoords,Area,Ktherm,source);
% check for any flux boundary conditions on the right side
%
Fbc = zeros(nNode,1);
if(elem==elemFlux)
Fbc(nodeFlux,1) = Area(Length)*flux;
end
% assemble this elements stiffness contribution into the global
%
for i=1:nNode;
for j=1:nNode;
row = connect(elem,i);
col = connect(elem,j);
K(row,col) = K(row,col) + Ke(i,j);
end
end
% assemble this elements force into the global force
%
for i=1:nNode;
row = connect(elem,i);
F(row) = F(row) + Fe(i) + Fbc(i);
end
end
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% impose temperature boundary condition using the penalty method
%
Ktmp = K(nodeT,nodeT); % needed to compute the reaction flux later
Ftmp = F(nodeT); % needed to compute the reaction flux later
%
% compute penalty parameter
%
penalty = (1.e6)*(trace(K)/nNodes);
%
% modify the stiffness and force
%
K(nodeT,nodeT) = penalty;
F(nodeT) = penalty*Tbc;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Global system solves
%
% compute the temperature field
%
T = K\F;
%
% compute the reaction flux
%
K(nodeT,nodeT) = Ktmp;
reaction = K(nodeT,:)*T - Ftmp;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Postprocessing to plot the temperature field
%
postprocess(Length,coord,T,reaction,nodeT)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
function postprocess(Length,coord,T,reaction,nodeT)
figure(1);
clf;
set(gca,'FontSize',12);
plot(coord,T,'ko-','LineWidth',2);
hold off;
xlim([0 Length]);
xlabel('Position');
ylabel('Temperature');
%
fprintf('Reaction flux at node %i is %d \n', nodeT, reaction);
% recall plotting the dT/dX one needs the B matrix in each element
end
function [Ke,Fe] = element(nInt,nNode,nodeCoords,Area,Ktherm,source)
% obtain gauss points and weights
%
if(nInt==1);
[xi,w] = GaussInt1Pt();
else
error('nInt is not programmed');
end
% obtain nodal coordinates
%
if(nNode==2);
x1 = nodeCoords(1,1);
x2 = nodeCoords(2,1);
else
error('nNode is not programmed');
end
% init
%
Ke = zeros(nNode,nNode);
Fe = zeros(nNode,1);
%
% loop over integration points
%
for intPt=1:nInt;
% compute shape functions and derivatives
%
if(nNode==2);
[N,B,Jac] = shapeLinear(x1,x2,xi(intPt));
else
error('nNode is not programmed');
end
% current location
%
x = N*nodeCoords;
% update the element stiffness
%
Ke = Ke + Jac*w(intPt)*transpose(B)*Area(x)*Ktherm(x)*B;
% update the element force vector
%
Fe = Fe + Jac*w(intPt)*transpose(N)*source(x);
end % loop over integration points
return;
end
function [xi,w] = GaussInt1Pt()
% Gauss integration locations and weights for 1pt integration
xi = 0.0;
%
w = 2.0;
return;
end
function [N,B,Jac] = shapeLinear(x1,x2,xi)
% shape functions and derivatives for a 2-node 1D element
% element length
%
Le = x2 - x1;
% the shape function matrix
%
N = (1.0/2.0)*[1.0-xi 1.0+xi];
% derivatives of shape functions
%
B = (1.0/Le)*[-1.0 1.0];
% the mapping jacobian
%
Jac = Le/2.0;
return;
end