**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 **

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-in^{4}. 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'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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-in^{4}. 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-in^{4}.

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 **

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.

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*