Computational Solid Mechanics Homework #2 Liu Yazhuo
3. FEM Analysis and coding (By MATLAB
®
R2018a - Academic use)
First, we have to do some preparation for FEM computing, such as element division, node area etc.
%% Area-computing function
linearArea = @ (x, A0, A1, L) (A1-A0)*x/L + A0;
%% Compute element information: location, average area etc.
n = m + 1; % Number of Nodes
x = linspace(0,L,n); % Compute node coordinates
elength = L/m; % element length
areaAtNode = linearArea(x, A0, A1, L);
areaAverage = ( areaAtNode(1:end-1) + areaAtNode(2:end) ) / 2;
nodeElementRelation = [1:m;1:m;2:n];
% element number in first row, the corresponding node number in
second and third row. every column is an element
Prepare stiffness matrix for each element and assemble them:
k_temp = [1,-1;-1,1];
for i = 1:m
k(:,:,i) = (E * areaAverage(i)/elength).*k_temp;
end
clear k_temp
totalStiffnessMatrix = AssemblyMatrix(k,nodeElementRelation,1);
AssemblyMatrix is a function to assemble the stiffness matrix by input a set of element stiffness
matrix and the corresponding element-node relationship. The function is as follows.
function totalStiffness = AssemblyMatrix(k,nodeElementRelation,n)
% input parameter format:
% k: 3D array and k(:,:,i) is the stiffness matrix for element i.
% nodeElementRelation is element number in first row,
% the corresponding node number in second and third row.
% n is the DOF of node
totalStiffness = zeros((size(nodeElementRelation,2)+1)*n);
for iCount = 1:size(nodeElementRelation,2)
einf = nodeElementRelation(:,iCount);
iNode = einf(2);
jNode = einf(3);