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

```% 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=0
%       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

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)
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-2*nu0)*[ 1-nu0      nu0        0;
nu0    1-nu0        0;
0        0  1/2-nu0 ];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GENERATE FINITE ELEMENT MESH
%
% Here we gnerate the finte element mesh (using the approriate elements).
% I won't go into too much detail about how to use these functions.  If
% one is interested one can type - help 'function name' at the matlab comand
% line to find out more about it.
%
% The folowing data structures are used to describe the finite element
% discretization:
%
%   node    - is a matrix of the node coordinates, i.e. node(I,j) -> x_Ij
%   element - is a matrix of element connectivities, i.e. the connectivity
%             of element e is given by > element(e,:) -> [n1 n2 n3 ...];
%
% To apply boundary conditions a description of the boundaries is needed.  To
% accomplish this we use a separate finite element discretization for each
% boundary.  For a 2D problem the boundary discretization is a set of 1D elements.
%
%   rightEdge -  a element connectivity matrix for the right edge
%   leftEdge  -  I'll give you three guesses
%
% These connectivity matricies refer to the node numbers defined in the
% coordinate matrix node.

disp([num2str(toc),'   GENERATING MESH'])
switch elemType
case 'Q4'           % here we generate the mesh of Q4 elements
nnx=numx+1;
nny=numy+1;
node=square_node_array([0 -c],[L -c],[L c],[0 c],nnx,nny);

inc_u=1;
inc_v=nnx;
node_pattern=[ 1 2 nnx+2 nnx+1 ];

element=make_elem(node_pattern,numx,numy,inc_u,inc_v);

case 'Q9'           % here we generate a mehs of Q9 elements
nnx=2*numx+1;
nny=2*numy+1;
node=square_node_array([0 -c],[L -c],[L c],[0 c],nnx,nny);

inc_u=2;
inc_v=2*nnx;
node_pattern=[ 1 3 2*nnx+3 2*nnx+1 2 nnx+3 2*nnx+2 nnx+1 nnx+2 ];

element=make_elem(node_pattern,numx,numy,inc_u,inc_v);

otherwise %'T3'    % and last but not least T3 elements
nnx=numx+1;
nny=numy+1;
node=square_node_array([0 -c],[L -c],[L c],[0 c],nnx,nny);

node_pattern1=[ 1 2 nnx+1 ];
node_pattern2=[ 2 nnx+2 nnx+1 ];
inc_u=1;
inc_v=nnx;

element=[make_elem(node_pattern1,numx,numy,inc_u,inc_v);
make_elem(node_pattern2,numx,numy,inc_u,inc_v) ];

end

% DEFINE BOUNDARIES
%   Here we define the boundary discretizations.
uln=nnx*(nny-1)+1;    % upper left node number
urn=nnx*nny;          % 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:2*nnx:(uln-1); (lrn+2*nnx):2*nnx:urn; (lrn+nnx):2*nnx:urn ]';
leftEdge =[ uln:-2*nnx:(lrn+1); (uln-2*nnx):-2*nnx: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 2*numnode. 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(2*numnode,1);          % nodal displacement vector
K=sparse(2*numnode,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

switch elemType  % define quadrature rule
case 'Q9'
otherwise
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

[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)/(2*I0);          % y traction at quadrature point
f(sctry)=f(sctry)+N*fyPt*detJ0*wt;  % scatter force into global force vector

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;

[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)/(2*I0);         % y traction at quadrature point
fxPt=P*L*yPt/I0;                    % x traction at quadrature point
f(sctry)=f(sctry)+N*fyPt*detJ0*wt;
f(sctrx)=f(sctrx)+N*fxPt*detJ0*wt;

end  % of element loop

% set the force at the nodes on the top and bottom edges to zero (traction free)
% TOP EDGE
topEdgeNodes = find(node(:,2)==c);  % finds nodes on the top edge
f(topEdgeNodes)=0;
f(topEdgeNodes+numnode)=0;

% BOTTOM EDGE
bottomEdgeNodes = find(node(:,2)==-c); % finds nodes on the bottom edge
f(bottomEdgeNodes)=0;
f(bottomEdgeNodes+numnode)=0;

%%%%%%%%%%%%%%%%%%%%% COMPUTE STIFFNESS MATRIX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp([num2str(toc),'   COMPUTING STIFFNESS MATRIX'])
switch elemType  % define quadrature rule
case 'Q9'
case 'Q4'
otherwise
end

for e=1:numelem                          % start of element loop

sctr=element(e,:);           % element scatter vector
sctrB=[ sctr sctr+numnode ]; % vector that scatters a B matrix

nn=length(sctr);

[N,dNdxi]=lagrange_basis(elemType,pt); % element shape functions

J0=node(sctr,:)'*dNdxi;                % element Jacobian matrix
invJ0=inv(J0);
dNdx=dNdxi*invJ0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% COMPUTE B MATRIX
%         _                                     _
%        |   N_1,x  N_2,x  ...     0      0  ... |
%  B  =  |       0      0  ... N_1,y  N_2,y  ... |
%        |   N_1,y  N_2,y  ... N_1,x  N_2,x  ... |
%         -                                     -
B=zeros(3,2*nn);
B(1,1:nn)       = dNdx(:,1)';
B(2,nn+1:2*nn)  = dNdx(:,2)';
B(3,1:nn)       = dNdx(:,2)';
B(3,nn+1:2*nn)  = dNdx(:,1)';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% COMPUTE ELEMENT STIFFNESS AT QUADRATURE POINT
K(sctrB,sctrB)=K(sctrB,sctrB)+B'*C*B*W(q)*det(J0);

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)=bcwt*speye(length(udofs)); % put ones*bcwt 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=dNdxi*invJ0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% COMPUTE B MATRIX
B=zeros(3,2*nn);
B(1,1:nn)       = dNdx(:,1)';
B(2,nn+1:2*nn)  = dNdx(:,2)';
B(3,1:nn)       = dNdx(:,2)';
B(3,nn+1:2*nn)  = dNdx(:,1)';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% COMPUTE ELEMENT STRAIN AND STRESS AT STRESS POINT
strain=B*U(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('************************************************')
```