Objectives

This lab exercise aims to perform a convegence study for 1-D deflection in a simply supported Euler-Bernoulli beam. For the first problem, the beam has a uniform cross section and experiences a uniformly applied load. Here a cubic Hermite finite element (FE) basis is used. For the second problem, the beam has a non-uniform cross section and experiences non-uniform loading conditions including point loads. Here a linear FE basis is used. Both problems are tested for asymptotic convergence rates in a convergence study that utilizes both the maximum displacement and energy norm convergence methods.

Problem Statement

PART A


Cubic Hermite FE Basis
Uniform Load, Uniform Cross-Section


Figure 1: 1-D deflection on a simply supported Euler-Bernoulli beam with uniform cross-section and uniform load.

The first lab problem involves a simply supported uniform beam, 120-in in length, under uniform loading as seen in Figure 1. The assumption of small displacements reduces the problem to one of 1-D, perfectly elastic deflections. The beam material is homogeneous with a modulus of elastisity of 30,000,000-psi. Also, the cross-section remains constant throughout the length of the beam resulting in a constant moment of intertia of 20.83333-in4. A load of 8.33333 lbf-in is applied uniformly across the entire length of the beam. Thus the beam is of Euler-Bernoulli type resulting in the following differential equation:

This fourth-order equation thus requires four boundary conditions representing the beam displacement and slope at L=0 and L=120-in. Since the beam is simply supported, all four conditions are of Dirichlet type and equation to zero. Also, the deflection variable, y, requires second order differentiation, which places a restriction upon the finite element (FE) basis choice. A cubic Hermite basis allows for second order differentiation and thus satisfies the restriction. With these developments in mind, the Matlab FEMLib template pseudo-code for the Galerkin Weak Statement becomes:

Finally, it is noted that the biharmonic derivative seen in the Euler-Bernoulli formulation requires a refinement of the convergence theory, namely:

The new variable, m, equals 2 for the fourth order derivative (it equalled 1 for the second order derivatives of the previous labs). Thus, for a cubic Hermite FE basis, k is equal to three resulting in an expected convergence slope of 4 for the problem solution.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                                        %
% Joseph Tipton                                                                          %
% ES551, LAB #6a                                                                         %
% 1-D, steady-state, deflection for  a simply supported Euler-Bernoulli beam problem     %   
% using the cubic Hermite element formulation                                            %
%                                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all;
clc;
global X1;

M = [2 4 8 16 32 64 128];                   % number of elements for each run

% set up physical parameters
E = 30e6;                                   % elastic modulus in psi
Izz = 20.83333;                             % moment of inertia in in^4
L = 120;                                    % length of beam in inches
p = 8.33333;                                % load in lbf/in

for i = 1:size(M,2)

    %set up mesh
    nnodes = M(i) + 1;
    X1 = linspace(0,L,nnodes);

	% set up nodal load data in Hermite form
	% need value and slope of load distribution at each node
	loadvals = ones(1,nnodes) * p;          % load values at each node
	slopevals = zeros(1,nnodes);            % slope of the load distribution at each node
	
	% loadvals and slopevals were formed as rows.  Concatenate them one on
	% top of the other, to form a 2 by nnodes array.
	PP = [loadvals;slopevals];
	% now stack the columns of PP on top of each other to form one long
	% column.  The first entry in P is the value of the load at node 1, the
	% second entry is the slope of the load at node 1, the third entry is the 
	% value of the load distribution at node 2, ect
	P = PP(:);
	
	load femHlib;                            % load cubic Hermite element library matrices
	
	K = ass1dH(E*Izz,[],[],-3,A222H,X1);     % form the global stiffness matrix
	b = ass1dH([],[],[],1,A200H,X1) * P;     % form the global load column matrix
	
	%  impose zero displacement boundary condtions at each end
	Kdiri = K;                               % copy K
	Kdiri(1,:) = zeros(1,2*nnodes);          % left end displacement
	Kdiri(1,1) = 1;                          % left end displacement
	b(1,1) = 0;                              % left end displacement
	Kdiri(2*nnodes-1,:) = zeros(1,2*nnodes); % right end displacement
	Kdiri(2*nnodes-1,2*nnodes-1) = 1;        % right end displacement
	b(2*nnodes-1,1) = 0;                     % right end displacement
	
	% no action required for zero moment at each end
	
	% solve the linear system 
	Y = Kdiri \ b;
    
    % save energy norm
    Enorm(i,1) = 1/2 * Y' * K * Y;
    
    % separate deflection and slope data
    for j = 1:size(X1,2)
        Deflection(j,1) = Y(2*j - 1);
        Moment(j,1) = Y(2*j);
    end
    
    % save maximum displacement value
    Ymax(i,1) = max(Deflection);
end

delta_Y = Ymax(2:end) - Ymax(1:end-1);
eh2 = delta_Y / (2^2 - 1);
Yexact = eh2 + Ymax(2:end);
slope_Taylor = log10(eh2(1:end-1) ./ eh2(2:end)) / log10(2);
error_est =  (Enorm(2:end) - Enorm(1:end-1)) ./  (2^2 - 1);
Enorm_exact_est = Enorm(2:end) + error_est;
slope_Enorm = log10(error_est(1:end-1) ./ error_est(2:end)) / log10(2);
format short				                 % put Matlab back to default format
answer(:,1) = M';
answer(:,2) = Ymax;
answer(:,3) = [0;delta_Y];
answer(:,4) = [0;eh2];
answer(:,5) = [0;Yexact];
answer(:,6) = [0;0;slope_Taylor];
answer(:,7) = Enorm;
answer(:,8) = [0;error_est];
answer(:,9) = [0;Enorm_exact_est];
answer(:,10) = [0;0;slope_Enorm];

fprintf('Mesh     Ymax  Ymax_change  TS_error  Ymax_exact  TS_Slope  Enorm ...
        ... Enorm_error  Enorm_exact  Enorm_Slope\n\n')
fprintf('%6.0f  %6.4f  %12.6g  %12.6g  %12.4f  %12.6g  %12.6g  %12.6g  %12.6g ...
        ... %12.6g\n',answer');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

PART B


Linear FE Basis
Non-Uniform Load, Non-Uniform Cross-Section


Figure 2: 1-D deflection on a simply supported Euler Bernoulli beam with non-uniform cross-section and non-uniform loading.

In the second lab problem, the Euler-Bernoulli beam is still simply supported, homogeneous with a constant modulus of elastisity of 30,000,000-psi, and still has zero displacement and slope Dirichlet BCs. Now, however, the cross-section and loading are non-uniform. As seen in Figure 2, a point load of 500-lbf is applied 30-in from the left end of the beam. Also, a linearly varying load is applied to the right side of the beam. At the midpoint, the load is zero and at the right end of the beam (L=120 in), the load is 100-lbf/in.


Figure 3: "T" beam cross section schematic.

Figure 3 illustrates how the beam cross-section varies. The entire beam consists of a 2x6-in rectangular beam which results in a moment of inertia of 4-in4. For the middle 1/4 of the beam, another 2x6-in rectangular beam is placed perpindicular as shown. Due to its orientation against the plane of deflection, this beam greatly increases the moment of inertia to 136-in4.

Thus the beam consists of two discrete moments of inertia applied as shown in Figure 2.

The method of solution for this problem involves a separation of the Euler-Bernoulli equation into two coupled second-order differential equations.

Now that the equations are reduced to second-order, any FE basis may be used and a Lagrange linear basis is chosen. The varying moment of inertia requires a hypermatrix on the deflection equation. The resulting MATLAB FEMLib pseudo-code of the Galerkin Weak Statment follows:

It should be noted that the point load must be placed at an element node. Thus an eight element mesh is required, at minimum, for a solution to proceed.

Finally, the non-smoothness of the beam loading and moment of inertia are expected to degrade the convergence performance of the solution. This is represented mathematically as follows:

where r, the order of analytical solution differentiability, is unknown.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                                        %
% Joseph Tipton                                                                          %
% ES551, LAB #6b                                                                         %
% 1-D, steady-state, deflection for  a simply supported Euler-Bernoulli beam problem     %   
% using a linear FE basis with non-uniform cross-section and non-uniform loading.        %
%                                                                                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all
clc
global X1;

Mesh = [8 16 32 64 128 256 512 1024];       % number of elements for each run

% set up physical parameters
E = 30e6;                                   % elastic modulus in psi
Izz_A = 4;                                  % "T" beam moment of inertia (in^4)
Izz_B = 134;                                % "I" beam moment of inertia (in^4)
L = 120;                                    % length of beam in inches
p = 100;                                    % load strength parameter (lbf/in)
P = 500;                                    % point load (lbf)

for i = 1:size(Mesh,2)

    %set up mesh
    nnodes = Mesh(i) + 1;
    X1 = linspace(0,L,nnodes);

    % set up nodal load data
    P_dist = (p * ((2/L)*X1 - 1))';
    P_dist(1:(Mesh(i)/2 + 1),1) = 0;
        
    load femlib;
	
    % set up M system matrix
    Msys =  asjac1D([],[],[],-1,A3011L,[]);
    Msysdiri = Msys;                        % copy Msys
    Msysdiri(1,:) = zeros(1,nnodes);        % apply diri to node 1
    Msysdiri(1,1) = 1;
    Msysdiri(nnodes,:) = zeros(1,nnodes);   % apply diri to last node
    Msysdiri(nnodes,nnodes) = 1;
	
    % set up M right hand side
    Mrhs =  asres1D([],[],[],1,A3000L,P_dist);
    Mrhs(1,1) = 0;                          % enforce diri condition at node 1
    Mrhs(nnodes,1) = 0;                     % enforce diri at last node
    point = 0.25*Mesh(i) + 1;               % locate node for point load
    Mrhs(point,1) = Mrhs(point,1) + P  ;    % enforce point load
	
    % solve for M
    M = Msysdiri \ Mrhs;
    Mmax(i,1) = max(M);                     % save Mmax
    node_Mmax(i,1) = find(M == max(M));     % save Mmax node location
    
    % set up nodal moment of inertia data
    Izz = ones(1,nnodes) * Izz_A;           % create nodally averaged Izz array
    left = (3 * Mesh(i) / 8) + 1;           % left-most node for plate addition
    right = (5 * Mesh(i) / 8) + 1;          % right-most node for plate addition
    Izz(1,left:right) = Izz_B;              % apply plate to middle 1/4 of beam
	
    % set up Y system matrix
    Ysys = asjac1D(E,[],Izz,-1,A3011L,[]);
    Ysysdiri = Ysys;
    Ysysdiri(1,:) = zeros(1,nnodes);        % apply diri to node 1
    Ysysdiri(1,1) = 1;
    Ysysdiri(nnodes,:) = zeros(1,nnodes);   % apply diri to last node
    Ysysdiri(nnodes,nnodes) = 1;
	
    % set up Y right hand side
    Yrhs = asres1D([],[],[],1,A3000L,M);
    Yrhs(1,1) = 0;                          % enforce diri condition at node 1
    Yrhs(nnodes,1) = 0;                     % enforce diri at last node
	
    % solve for Y
    Y = Ysysdiri \ Yrhs;
    Ymax(i,1) = max(Y);                     % save Ymax
    node_Ymax(i,1) = find(Y == max(Y));     % save Ymax node location
	
    % calculate the energy norms
    format long
    enormM(i,1) = 1/2 * M' * Msys * M;
    enormY(i,1) = 1/2 * Y' * Ysys * Y;
    
end

delta_Ymax = Ymax(2:end) - Ymax(1:end-1);
slope_Ymax = log10(delta_Ymax(1:end-1) ./ delta_Ymax(2:end)) / log10(2);
delta_Mmax = Mmax(2:end) - Mmax(1:end-1);
slope_Mmax = log10(delta_Mmax(1:end-1) ./ delta_Mmax(2:end)) / log10(2);

MAX(:,1) = Mesh';
MAX(:,2) = Mmax;
MAX(:,3) = node_Mmax;
MAX(:,4) = [0;delta_Mmax];
MAX(:,5) = [0;0;slope_Mmax];
MAX(:,6) = Ymax;
MAX(:,7) = node_Ymax;
MAX(:,8) = [0;delta_Ymax];
MAX(:,9) = [0;0;slope_Ymax];

fprintf('Mesh     Mmax  (Mmax)node  (Mmax)change  (Mmax)slope     Ymax  (Ymax)node ...
        ... (Ymax)change  (Ymax)slope\n\n')
fprintf('%6.0f  %6.4f  %12.6g  %12.6g  %12.6g  %6.4f  %12.6g  %12.6g  %12.6g\n',MAX');

delta_enormM =  enormM(2:end) - enormM(1:end-1);
slope_enormM = log10(delta_enormM(1:end-1) ./ delta_enormM(2:end)) / log10(2);
delta_enormY =  enormY(2:end) - enormY(1:end-1);
slope_enormY = log10(delta_enormY(1:end-1) ./ delta_enormY(2:end)) / log10(2);

ENORM(:,1) = Mesh';
ENORM(:,2) = enormM;
ENORM(:,3) = [0;delta_enormM];
ENORM(:,4) = [0;0;slope_enormM];
ENORM(:,5) = enormM;
ENORM(:,6) = [0;delta_enormM];
ENORM(:,7) = [0;0;slope_enormM];

fprintf('\n\nMesh        ||M||      change||M||     Slope(M)       ||Y|| ...
        ... change||Y||  Slope(Y)\n\n')
fprintf('%6.0f  %12.6g  %12.6g  %12.6g  %12.6g  %12.6g  %12.6g\n',ENORM');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Discussion of Results

PART A


Cubic Hermite FE Basis
Uniform Load, Uniform Cross-Section


Figure 1: Steady-state deflections across a simply supported Euler-Bernoulli beam with uniform load and uniform cross-section. Solved with 128 mesh elements with a cubic Hermite FE basis.

 


Figure 2: Steady-state moment distribution across a simply supported Euler-Bernoulli beam with uniform load and uniform cross-section. Solved with 128 mesh elements with a cubic Hermite FE basis.

Figures 1 and 2 give the displacement and slope solutions for the beam. As expected, Figure 1 shows the maximum displacement occurs at the center of the beam with the Dirichlet conditions of zero deflection maintained at either end. The maximum displacement is approximatly 0.036-in which validates the assumption of small beam deflections.

Table 1: Taylor Series truncation error convergence results for a cubic Hermite Finite Element basis.
=====================================================================
  Mesh    Ymax    Ymax_change     Truncation  Ymax_exact    TS_Slope 
  Size				  Error Est.      (Est.)             
_____________________________________________________________________
     2  0.0360              0              0      0.0000           0 
     4  0.0360  -1.87350e-016  -6.24500e-017      0.0360           0 
     8  0.0360  -4.85723e-016  -1.61908e-016      0.0360     -1.3744 
    16  0.0360   7.57727e-015   2.52576e-015      0.0360    -3.96347 
    32  0.0360   3.20438e-013   1.06813e-013      0.0360    -5.40222 
    64  0.0360  -4.00534e-013  -1.33511e-013      0.0360   -0.321878 
   128  0.0360   5.35544e-014   1.78515e-014      0.0360     2.90285 
=====================================================================

Table 1 gives the convergence results of the mesh refinement based upon the Taylor Series truncation error at the maximum beam displacement. As expected, the simple case of a simply supported beam with uniform cross section and uniform loading results in a nodally exact solution. The slope calculations then become useless due to computer round-off effects as seen.

Table 2: Energy Norm convergence results for a cubic Hermite Finite Element basis.
==========================================================
  Mesh    Enorm   Enorm_error   Enorm_exact   Enorm_Slope 
__________________________________________________________
     2  11.3999             0             0             0 
     4  11.5124     0.0374998       11.5499             0 
     8  11.5195    0.00234374       11.5218             4 
    16  11.5199   0.000146483         11.52             4 
    32  11.5199  9.15528e-006       11.5199       3.99999 
    64  11.5199  5.72081e-007       11.5199       4.00031 
   128  11.5199    3.659e-008       11.5199        3.9667 
==========================================================

Table 2 gives the convergence results of the mesh refinement based upon the energy norm. Even though the solutions are nodally exact, the energy norm method retains the ability to resolve the error. The resulting slope calculations indicate a value of 4 that agrees well with the expected theory.

PART B


Linear FE Basis
Non-Uniform Load, Non-Uniform Cross-Section


Figure 3: Steady-state deflections across a simply supported Euler-Bernoulli beam with non-uniform load and non-uniform cross-section. Solved with 1024 mesh linear basis elements.

 


Figure 4: Steady-state moment distribution across a simply supported Euler-Bernoulli beam with non-uniform load and non-uniform cross-section. Solved with 1024 mesh linear basis elements.

Figures 3 and 4 give the displacement and slope solutions for the beam. Figure 3 verifies the assumption of small beam deflections and shows the zero displacement Dirichlet conditiosn at either end of the beam. The effect of the much higher moment of inertia applied at the middle quarter of the beam is seen to constrain any deflection and results in the "cap" seen in the figure. Figure 4 shows zero moments at either end of the beam in keeping with the specified Dirichlet BCs. Also, the effect of the point load, applied at L=30-in, is seen by an abrupt change in the slope of the moment.

Table 3: Maximum moment and displacement convergence results for a linear FE basis.
=============================================================================================
  Mesh        Mmax   Mmax           Mmax      Mmax      Ymax   Ymax          Ymax      Ymax  
	  (lbf-in)   node         change     slope  (lbf-in)   node        change     slope  
_____________________________________________________________________________________________
     8  42187.5000      6              0         0    0.3718      5             0         0  
    16  42773.4375     12        585.937         0    0.4320     10     0.0601655         0  
    32  42773.4375     23  -1.45519e-011   45.1946    0.4554     18     0.0234014   1.36234  
    64  42797.2412     44        23.8037  -40.5731    0.4654     36    0.00999951   1.22666  
   128  42801.1322     88        3.89099   2.61298    0.4699     71    0.00455891   1.13317  
   256  42803.0491    174        1.91689   1.02137    0.4721    140    0.00217137   1.07008  
   512  42803.0670    348      0.0178814   6.74416    0.4731    279    0.00105805    1.0372  
  1024  42803.3008    694       0.233799  -3.70874    0.4737    557   0.000522196   1.01875  
=============================================================================================

Table 3 gives the convergence results of the mesh refinement in terms of the maximum displacement. Since the data non-smoothness introduces an unknown into the convergence theory, the error can no longer be approximated. Instead, the change in maximum displacement is used as a measure of accuracy to determine when the solution falls within 0.1%. As seen, this occurs at a mesh size of 1,024 elements. The convergence slope for the moment solution oscillates and the cause is unknown. The displacement convergence slope, however, clearly indicates a linear convergence rate.

Table 3: Energy Norm convergence results for a linear FE basis.
=========================================================================================
  Mesh         ||M||         ||M||     ||M||          ||Y||         ||Y||         ||Y||  
			    change     slope	                   change         slope  
_________________________________________________________________________________________
     8  4.74609e+007             0         0   4.74609e+007             0             0  
    16  4.88452e+007  1.38428e+006         0   4.88452e+007  1.38428e+006             0  
    32  4.91954e+007        350189   1.98293   4.91954e+007        350189       1.98293  
    64  4.92832e+007       87804.8   1.99576   4.92832e+007       87804.8       1.99576  
   128  4.93052e+007       21967.3   1.99894   4.93052e+007       21967.3       1.99894  
   256  4.93107e+007       5492.83   1.99974   4.93107e+007       5492.83       1.99974  
   512   4.9312e+007       1373.27   1.99993    4.9312e+007       1373.27       1.99993  
  1024  4.93124e+007       343.321   1.99998   4.93124e+007       343.321       1.99998  
=========================================================================================

Finally, Table 4 lists the convergence results of the mesh refinement in terms of the energy norm. Both the moment and displacement solutions indicate a convergence slope of 2. This is puzzling since the theory would indicate a maximum slope of 2 for a linear FE basis with further degredation expected due to data non-smoothness. The author is at a loss to explain this occurrance.

Conclusion

In conclusion, a convergence study has been completed for a simply supported Euler-Bernoulli beam. In the first part, the beam is homogeneous and uniform with a uniform applied load. This results in nodally exact solutions. A cubic Hermite FE basis yields convergence slopes in the energy norm method near 4. For the second part, the beam is homogeneous with a non-uniform cross-section and non-uniform loading. The resulting data non-smoothless results in a degredation of the convergence rate. Reformulating the fourth-order equation of deflection into two coupled second order differential equations (in moment and deflection) allows the use of Lagrange linear FE basis.

Grade

10.18.04 Joe, your lab 6 is done very professionally, throughout. Results stand to reason except for the slope of 2 for the non-uniform loading case. How you got two must be associated with the I distribution you chose?? Most importantly is the issue of convergence. If the results converge then you can estimate the mesh necessary for engineering accuracy independent of the convergence rate. Conclusions are right on for your data. Overall excellent job!! Score is 7/6 AJB

ML>