Download Programing the Finite Element Method with Matlab and more Study Guides, Projects, Research Computer Numerical Control in PDF only on Docsity!
Programing the Finite Element Method with
Matlab
Jack Chessa∗
3rd October 2002
1 Introduction
The goal of this document is to give a very brief overview and direction in the writing of finite element code using Matlab. It is assumed that the reader has a basic familiarity with the theory of the finite element method, and our attention will be mostly on the implementation. An example finite element code for analyzing static linear elastic problems written in Matlab is presented to illustrate how to program the finite element method. The example program and supporting files are available at http://www.tam.northwestern.edu/jfc795/Matlab/
1.1 Notation
For clarity we adopt the following notation in this paper; the bold italics font v denotes a vector quantity of dimension equal to the spacial dimension of the problem i.e. the displacement or velocity at a point, the bold non-italicized font d denotes a vector or matrix which is of dimension of the number of unknowns in the discrete system i.e. a system matrix like the stiffness matrix, an uppercase subscript denotes a node number whereas a lowercase subscript in general denotes a vector component along a Cartesian unit vector. So, if d is the system vector of nodal unknowns, uI is a displacement vector of node I and uIi is the component of the displacement at node I in the i direction, or uI · ei. Often Matlab syntax will be intermixed with mathematical notation
∗Graduate Research Assistant, Northwestern University ([email protected])
which hopefully adds clarity to the explanation. The typewriter font, font, is used to indicate that Matlab syntax is being employed.
2 A Few Words on Writing Matlab Programs
The Matlab programming language is useful in illustrating how to program the finite element method due to the fact it allows one to very quickly code numerical methods and has a vast predefined mathematical library. This is also due to the fact that matrix (sparse and dense), vector and many linear algebra tools are already defined and the developer can focus entirely on the implementation of the algorithm not defining these data structures. The extensive mathematics and graphics functions further free the developer from the drudgery of developing these functions themselves or finding equivalent pre-existing libraries. A simple two dimensional finite element program in Matlab need only be a few hundred lines of code whereas in Fortran or C++ one might need a few thousand. Although the Matlab programming language is very complete with re- spect to it’s mathematical functions there are a few finite element specific tasks that are helpful to develop as separate functions. These have been programed and are available at the previously mentioned web site. As usual there is a trade off to this ease of development. Since Matlab is an interpretive language; each line of code is interpreted by the Matlab command line interpreter and executed sequentially at run time, the run times can be much greater than that of compiled programming languages like Fortran or C++. It should be noted that the built-in Matlab functions are already compiled and are extremely efficient and should be used as much as possible. Keeping this slow down due to the interpretive nature of Matlab in mind, one programming construct that should be avoided at all costs is the for loop, especially nested for loops since these can make a Matlab programs run time orders of magnitude longer than may be needed. Often for loops can be eliminated using Matlab’s vectorized addressing. For example, the following Matlab code which sets the row and column of a matrix A to zero and puts one on the diagonal
for i=1:size(A,2) A(n,i)=0; end for i=1:size(A,1) A(i,n)=0; end
4 Finite Element Data Structures in Matlab
Here we discuss the data structures used in the finite element method and specifically those that are implemented in the example code. These are some- what arbitrary in that one can imagine numerous ways to store the data for a finite element program, but we attempt to use structures that are the most flexible and conducive to Matlab. The design of these data structures may be depend on the programming language used, but usually are not significantly different than those outlined here.
4.1 Nodal Coordinate Matrix
Since we are programming the finite element method it is not unexpected that we need some way of representing the element discretization of the domain. To do so we define a set of nodes and a set of elements that connect these nodes in some way. The node coordinates are stored in the nodal coordinate matrix. This is simply a matrix of the nodal coordinates (imagine that). The dimension of this matrix is nn × sdim where nn is the number of nodes and sdim is the number of spacial dimensions of the problem. So, if we consider a nodal coordinate matrix nodes the y-coordinate of the nth^ node is nodes(n,2). Figure 1 shows a simple finite element discretization. For this simple mesh the nodal coordinate matrix would be as follows:
nodes =
4.2 Element Connectivity Matrix
The element definitions are stored in the element connectivity matrix. This is a matrix of node numbers where each row of the matrix contains the con- nectivity of an element. So if we consider the connectivity matrix elements that describes a mesh of 4-node quadrilaterals the 36th element is defined by the connectivity vector elements(36,:) which for example may be [ 36 42 13 14] or that the elements connects nodes 36 → 42 → 13 → 14. So for
the simple mesh in Figure 1 the element connectivity matrix is
elements =
.^ (2)
Note that the elements connectivities are all ordered in a counter-clockwise fashion; if this is not done so some Jacobian’s will be negative and thus can cause the stiffnesses matrix to be singular (and obviously wrong!!!).
4.3 Definition of Boundaries
In the finite element method boundary conditions are used to either form force vectors (natural or Neumann boundary conditions) or to specify the value of the unknown field on a boundary (essential or Dirichlet boundary conditions). In either case a definition of the boundary is needed. The most versatile way of accomplishing this is to keep a finite element discretization of the necessary boundaries. The dimension of this mesh will be one order less that the spacial dimension of the problem (i.e. a 2D boundary mesh for a 3D problem, 1D boundary mesh for a 2D problem etc. ). Once again let’s consider the simple mesh in Figure 1. Suppose we wish to apply a boundary condition on the right edge of the mesh then the boundary mesh would be the defined by the following element connectivity matrix of 2-node line elements
right Edge =
[
]
Note that the numbers in the boundary connectivity matrix refer to the same node coordinate matrix as do the numbers in the connectivity matrix of the interior elements. If we wish to apply an essential boundary conditions on this edge we need a list of the node numbers on the edge. This can be easily done in Matlab with the unique function. nodesOnBoundary = unique(rightEdge); This will set the vector nodesOnBoundary equal to [2 4 6]. If we wish to from a force vector from a natural boundary condition on this edge we simply loop over the elements and integrate the force on the edge just as we would integrate any finite element operators on the domain interior i.e. the stiffness matrix K.
where nn is again the number of nodes in the discretization. This mapping is n = sdim(I − 1) + i. (10)
With this mapping the i component of the displacement at node I is located as follows in d dIi = d( sdim*(I-1) + i ). (11) The other option is to group all the like components of the nodal un- knowns in a contiguous portion of d or as follows
d =
u 1 x u 2 x .. . unx u 1 y u 2 y .. .
The mapping in this case is
n = (i − 1)nn + I (13)
So for this structure the i component of the displacement at node I is located at in d dIi = d( (i-1)*nn + I ) (14)
For reasons that will be appreciated when we discuss the scattering of element operators into system operators we will adopt the latter dof mapping. It is important to be comfortable with these mappings since it is an operation that is performed regularly in any finite element code. Of course which ever mapping is chosen the stiffness matrix and force vectors should have the same structure.
5 Computation of Finite Element Operators
At the heart of the finite element program is the computation of finite element operators. For example in a linear static code they would be the stiffness matrix
K =
Ω
BT^ C B dΩ (15)
and the external force vector
f ext^ =
Γt
Nt dΓ. (16)
The global operators are evaluated by looping over the elements in the dis- cretization, integrating the operator over the element and then to scatter the local element operator into the global operator. This procedure is written mathematically with the Assembly operator A
K = Ae
Ωe
BeT^ C Be^ dΩ (17)
5.1 Quadrature
The integration of an element operator is performed with an appropriate quadrature rule which depends on the element and the function being inte- grated. In general a quadrature rule is as follows
∫ (^) ξ=
ξ=− 1
f (ξ)dξ =
q
f (ξq)Wq (18)
where f (ξ) is the function to be integrated, ξq are the quadrature points and Wq the quadrature weights. The function quadrature generates a vector of quadrature points and a vector of quadrature weights for a quadrature rule. The syntax of this function is as follows
[quadWeights,quadPoints] = quadrature(integrationOrder, elementType,dimensionOfQuadrature);
so an example quadrature loop to integrate the function f = x^3 on a trian- gular element would be as follows
[qPt,qWt]=quadrature(3,’TRIANGULAR’,2); for q=1:length(qWt) xi = qPt(q); % quadrature point % get the global coordinte x at the quadrature point xi % and the Jacobian at the quadrature point, jac ... f_int = f_int + x^3 * jac*qWt(q); end
this involves modifying the system
Kd = f (19)
so that the essential boundary condition
dn = d¯n (20)
is satisfied while retaining the original finite element equations on the un- constrained dofs. In (20) the subscript n refers to the index of the vector d not to a node number. An easy way to enforce (20) would be to modify nth row of the K matrix so that
Knm = δnm ∀m ∈ { 1 , 2... N } (21)
where N is the dimension of K and setting
fn = d¯n. (22)
This reduces the nth^ equation of (19) to (20). Unfortunately, this destroys the symmetry of K which is a very important property for many efficient linear solvers. By modifying the nth^ column of K as follows
Km,n = δnm ∀m ∈ { 1 , 2... N }. (23)
We can make the system symmetric. Of course this will modify every equa- tion in (19) unless we modify the force vector f
fm = Kmn d¯n. (24)
If we write the modified kth^ equation in (19)
Kk 1 d 1 + Kk 2 d 2 +... Kk(n−1)dn− 1 + Kk(n+1)dn+1 +... + KkN dN = fk − Kkn d¯n (25)
it can be seen that we have the same linear equations as in (19), but just with the internal force from the constrained dof. This procedure in Matlab i s as follows
f = f - K(:,fixedDofs)fixedDofValues; K(:,fixedDofs) = 0; K(fixedDofs,:) = 0; K(fixedDofs,fixedDofs) = bcwtspeye(length(fixedDofs)); f(fixedDofs) = bcwt*fixedDofValues;
where fixedDofs is a vector of the indicies in d that are fixed, fixedDofValues is a vector of the values that fixedDofs are assigned to and bcwt is a weighing factor to retain the conditioning of the stiffness matrix (typically bcwt = trace(K)/N ).
6 Where To Go Next
Hopefully this extremely brief overview of programming simple finite element methods with Matlab has helped bridge the gap between reading the theory of the finite element method and sitting down and writing ones own finite element code. The examples in the Appendix should be looked at and run, but also I would suggest trying to write a simple 1D or 2D finite element code from scratch to really solidify the method in ones head. The examples can then be used as a reference to diminish the struggle. Good Luck!
B Example: Beam Bending Problem
The first example program solves the static bending of a linear elastic beam. The configuration of the problem is shown in Figure 3 and the program can be found at http://www.tam.northwestern.edu/jfc795/Matlab/ Examples/Static/beam.m The exact solution for this problem is as follows
σ 11 = −
P (L − x)y I σ 22 = 0 σ 12 =
P
2 I
(c^2 − y^2 )
u 1 = −
P y 6 EI
L^2 − (L − x)^2
u 2 =
P y 6 EI
(L − x)^3 − L^3
[
(4 + 5ν)c^2 + 3L^2
]
x + 3ν(L − x)y^2
This problem can be run with three element types; three node triangle element, a four node quadrilateral element and a nine node quadrilateral ele- ment. Also, one can choose between plane strain or plane stress assumption.
% beam.m % % Solves a linear elastic 2D beam problem ( plane stress or strain ) % with several element types. % % ^ y % | % --------------------------------------------- % | | % | | % ---------> x | 2c % | | % | L | % --------------------------------------------- % % with the boundary following conditions: % % u_x = 0 at (0,0), (0,-c) and (0,c) % u_y = 0 at (0,0) % % t_x = y along the edge x= % t_y = P*(x^2-c^2) along the edge x=L % % ****************************************************************************** % % This file and the supporting matlab files can be found at % http://www.tam.northwestern.edu/jfc795/Matlab % % by Jack Chessa % Northwestern University
% % ******************************************************************************
clear colordef black state = 0;
% ****************************************************************************** % *** I N P U T *** % ****************************************************************************** tic; disp(’************************************************’) disp(’*** S T A R T I N G R U N ***’) disp(’************************************************’)
disp([num2str(toc),’ START’])
% MATERIAL PROPERTIES E0 = 10e7; % Young’s modulus nu0 = 0.30; % Poisson’s ratio
% BEAM PROPERTIES L = 16; % length of the beam c = 2; % the distance of the outer fiber of the beam from the mid-line
% MESH PROPERTIES elemType = ’Q9’; % the element type used in the FEM simulation; ’T3’ is for a % three node constant strain triangular element, ’Q4’ is for % a four node quadrilateral element, and ’Q9’ is for a nine % node quadrilateral element.
numy = 4; % the number of elements in the x-direction (beam length) numx = 18; % and in the y-direciton. plotMesh = 1; % A flag that if set to 1 plots the initial mesh (to make sure % that the mesh is correct) % TIP LOAD P = -1; % the peak magnitude of the traction at the right edge
% STRESS ASSUMPTION stressState=’PLANE_STRESS’; % set to either ’PLANE_STRAIN’ or "PLANE_STRESS’ % nuff said.
% ****************************************************************************** % *** P R E - P R O C E S S I N G *** % ****************************************************************************** I0=2*c^3/3; % the second polar moment of inertia of the beam cross-section.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % COMPUTE ELASTICITY MATRIX if ( strcmp(stressState,’PLANE_STRESS’) ) % Plane Strain case C=E0/(1-nu0^2)[ 1 nu0 0; nu0 1 0; 0 0 (1-nu0)/2 ]; else % Plane Strain case C=E0/(1+nu0)/(1-2nu0)*[ 1-nu0 nu0 0; nu0 1-nu0 0; 0 0 1/2-nu0 ]; end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % GENERATE FINITE ELEMENT MESH %
% Here we define the boundary discretizations. uln=nnx(nny-1)+1; % upper left node number urn=nnxnny; % upper right node number lrn=nnx; % lower right node number lln=1; % lower left node number cln=nnx(nny-1)/2+1; % node number at (0,0) switch elemType case ’Q9’ rightEdge=[ lrn:2nnx:(uln-1); (lrn+2nnx):2nnx:urn; (lrn+nnx):2nnx:urn ]’; leftEdge =[ uln:-2nnx:(lrn+1); (uln-2nnx):-2nnx:1; (uln-nnx):-2*nnx:1 ]’; edgeElemType=’L3’;
otherwise % same discretizations for Q4 and T3 meshes rightEdge=[ lrn:nnx:(uln-1); (lrn+nnx):nnx:urn ]’; leftEdge =[ uln:-nnx:(lrn+1); (uln-nnx):-nnx:1 ]’; edgeElemType=’L2’;
end
% GET NODES ON DISPLACEMENT BOUNDARY % Here we get the nodes on the essential boundaries fixedNodeX=[uln lln cln]’; % a vector of the node numbers which are fixed in % the x direction fixedNodeY=[cln]’; % a vector of node numbers which are fixed in % the y-direction
uFixed=zeros(size(fixedNodeX)); % a vector of the x-displacement for the nodes % in fixedNodeX ( in this case just zeros ) vFixed=zeros(size(fixedNodeY)); % and the y-displacements for fixedNodeY
numnode=size(node,1); % number of nodes numelem=size(element,1); % number of elements
% PLOT MESH if ( plotMesh ) % if plotMesh==1 we will plot the mesh clf plot_mesh(node,element,elemType,’g.-’); hold on plot_mesh(node,rightEdge,edgeElemType,’bo-’); plot_mesh(node,leftEdge,edgeElemType,’bo-’); plot(node(fixedNodeX,1),node(fixedNodeX,2),’r>’); plot(node(fixedNodeY,1),node(fixedNodeY,2),’r^’); axis off axis([0 L -c c]) disp(’(paused)’) pause end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % DEFINE SYSTEM DATA STRUCTURES % % Here we define the system data structures % U - is vector of the nodal displacements it is of length 2numnode. The % displacements in the x-direction are in the top half of U and the % y-displacements are in the lower half of U, for example the displacement % in the y-direction for node number I is at U(I+numnode) % f - is the nodal force vector. It’s structure is the same as U, % i.e. f(I+numnode) is the force in the y direction at node I % K - is the global stiffness matrix and is structured the same as with U and f % so that K_IiJj is at K(I+(i-1)numnode,J+(j-1)numnode) disp([num2str(toc),’ INITIALIZING DATA STRUCTURES’]) U=zeros(2numnode,1); % nodal displacement vector
f=zeros(2numnode,1); % external load vector K=sparse(2numnode,2*numnode); % stiffness matrix
% a vector of indicies that quickly address the x and y portions of the data % strtuctures so U(xs) returns U_x the nodal x-displacements xs=1:numnode; % x portion of u and v vectors ys=(numnode+1):2*numnode; % y portion of u and v vectors
% ****************************************************************************** % *** P R O C E S S I N G *** % ******************************************************************************
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % COMPUTE EXTERNAL FORCES % integrate the tractions on the left and right edges disp([num2str(toc),’ COMPUTING EXTERNAL LOADS’])
switch elemType % define quadrature rule case ’Q9’ [W,Q]=quadrature( 4, ’GAUSS’, 1 ); % four point quadrature otherwise [W,Q]=quadrature( 3, ’GAUSS’, 1 ); % three point quadrature end
% RIGHT EDGE for e=1:size(rightEdge,1) % loop over the elements in the right edge
sctr=rightEdge(e,:); % scatter vector for the element sctrx=sctr; % x scatter vector sctry=sctrx+numnode; % y scatter vector for q=1:size(W,1) % quadrature loop pt=Q(q,:); % quadrature point wt=W(q); % quadrature weight [N,dNdxi]=lagrange_basis(edgeElemType,pt); % element shape functions J0=dNdxi’node(sctr,:); % element Jacobian detJ0=norm(J0); % determiniat of jacobian yPt=N’node(sctr,2); % y coordinate at quadrature point fyPt=P(c^2-yPt^2)/(2I0); % y traction at quadrature point f(sctry)=f(sctry)+NfyPtdetJ0*wt; % scatter force into global force vector end % of quadrature loop end % of element loop
% LEFT EDGE for e=1:size(leftEdge,1) % loop over the elements in the left edge
sctr=rightEdge(e,:); sctrx=sctr; sctry=sctrx+numnode; for q=1:size(W,1) % quadrature loop pt=Q(q,:); % quadrature point wt=W(q); % quadrature weight [N,dNdxi]=lagrange_basis(edgeElemType,pt); % element shape functions J0=dNdxi’node(sctr,:); % element Jacobian detJ0=norm(J0); % determiniat of jacobian yPt=N’node(sctr,2); fyPt=-P(c^2-yPt^2)/(2I0); % y traction at quadrature point fxPt=PLyPt/I0; % x traction at quadrature point
end % of element loop %%%%%%%%%%%%%%%%%%% END OF STIFFNESS MATRIX COMPUTATION %%%%%%%%%%%%%%%%%%%%%%
% APPLY ESSENTIAL BOUNDARY CONDITIONS disp([num2str(toc),’ APPLYING BOUNDARY CONDITIONS’]) bcwt=mean(diag(K)); % a measure of the average size of an element in K % used to keep the conditioning of the K matrix udofs=fixedNodeX; % global indecies of the fixed x displacements vdofs=fixedNodeY+numnode; % global indecies of the fixed y displacements
f=f-K(:,udofs)uFixed; % modify the force vector f=f-K(:,vdofs)vFixed; f(udofs)=uFixed; f(vdofs)=vFixed;
K(udofs,:)=0; % zero out the rows and columns of the K matrix K(vdofs,:)=0; K(:,udofs)=0; K(:,vdofs)=0; K(udofs,udofs)=bcwtspeye(length(udofs)); % put onesbcwt on the diagonal K(vdofs,vdofs)=bcwt*speye(length(vdofs));
% SOLVE SYSTEM disp([num2str(toc),’ SOLVING SYSTEM’]) U=K\f;
%****************************************************************************** %*** P O S T - P R O C E S S I N G *** %****************************************************************************** % % Here we plot the stresses and displacements of the solution. As with the % mesh generation section we don’t go into too much detail - use help % ’function name’ to get more details.
disp([num2str(toc),’ POST-PROCESSING’])
dispNorm=L/max(sqrt(U(xs).^2+U(ys).^2)); scaleFact=0.1*dispNorm; fn=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % PLOT DEFORMED DISPLACEMENT PLOT figure(fn) clf plot_field(node+scaleFact[U(xs) U(ys)],element,elemType,U(ys)); hold on plot_mesh(node+scaleFact[U(xs) U(ys)],element,elemType,’g.-’); plot_mesh(node,element,elemType,’w--’); colorbar fn=fn+1; title(’DEFORMED DISPLACEMENT IN Y-DIRECTION’)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % COMPUTE STRESS stress=zeros(numelem,size(element,2),3);
switch elemType % define quadrature rule case ’Q9’ stressPoints=[-1 -1;1 -1;1 1;-1 1;0 -1;1 0;0 1;-1 0;0 0 ]; case ’Q4’ stressPoints=[-1 -1;1 -1;1 1;-1 1]; otherwise
stressPoints=[0 0;1 0;0 1]; end
for e=1:numelem % start of element loop
sctr=element(e,:); sctrB=[sctr sctr+numnode]; nn=length(sctr); for q=1:nn pt=stressPoints(q,:); % stress point [N,dNdxi]=lagrange_basis(elemType,pt); % element shape functions J0=node(sctr,:)’dNdxi; % element Jacobian matrix invJ0=inv(J0); dNdx=dNdxiinvJ0; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % COMPUTE B MATRIX B=zeros(3,2nn); B(1,1:nn) = dNdx(:,1)’; B(2,nn+1:2nn) = dNdx(:,2)’; B(3,1:nn) = dNdx(:,2)’; B(3,nn+1:2nn) = dNdx(:,1)’; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % COMPUTE ELEMENT STRAIN AND STRESS AT STRESS POINT strain=BU(sctrB); stress(e,q,:)=C*strain; end end % of element loop
stressComp=1; figure(fn) clf plot_field(node+scaleFact[U(xs) U(ys)],element,elemType,stress(:,:,stressComp)); hold on plot_mesh(node+scaleFact[U(xs) U(ys)],element,elemType,’g.-’); plot_mesh(node,element,elemType,’w--’); colorbar fn=fn+1; title(’DEFORMED STRESS PLOT, BENDING COMPONENT’) %print(fn,’-djpeg90’,[’beam_’,elemType,’_sigma’,num2str(stressComp),’.jpg’])
disp([num2str(toc),’ RUN FINISHED’]) % *************************************************************************** % *** E N D O F P R O G R A M *** % *************************************************************************** disp(’************************************************’) disp(’*** E N D O F R U N ***’) disp(’************************************************’)
