function [] = Two_span_three_L3_C1_beam_react () % revised 3/9 % Requires: addpath /clear/www/htdocs/mech517/Akin_FEA_Lib % EI=1 w=-2,400 lb/ft P=-10,000 lb % (1) (2) | (3) % wwwwwwwwwwwww roller V % roller *-----*-----*--------*--------*---*---* fixed % NODE 1 2 3 4 5 6 7 % ^ ^ ^ % R 9,957 | | 17,226 lb 6,816 | | % -23,121 <--/ ft lb % X 0 5 10 16 22 25 28 ft. % % DOF 1,2 3,4 5,6 7,8 9,10 11,12 13,14 % EBC_Flag 10 00 10 00 00 00 11 % EBC_value 0. 0. 0.,0. fprintf('EI=1 w=-2,400 lb/ft P=-10,000 lb \n') fprintf('lb/ft^2 (1) (2) | (3) \n') fprintf(' wwwwwwwwwwwww V \n') fprintf('roller *-----*-----*--------*------*---*---* fixed\n') fprintf('NODE 1 2 3,roller 4 5 6 7 \n') fprintf('X 0 5 10 16 22 25 28 ft \n') fprintf('DOF 1,2 3,4 5,6 7,8 9,10 11,12 13,14 \n') fprintf('Reaction results are: \n') fprintf('R1=9,957 R3=17,226 R7=6,816 lb. M7=-23,121 ft lb \n') % problem controls n_g = 2 ; % number of DOF per node n_e = 3 ; % number of elements n_m = 7 ; % number of nodes in the mesh n_n = 3 ; % number of nodes per element n_d = n_g * n_m ; % number of system equations n_i = n_g * n_n ; % nmber of element DOF % allocate input arrays x = zeros (n_m, 1) ; nodes = zeros (n_e, n_n) ; % inputs nodes = [1 2 3 ; % element 1 connections 3 4 5 ; % element 2 connections 5 6 7 ] ; % element 3 connections x = [0 5 10 16 22 25 28]' ; % system x-coordinates EIs = [1 1 1] ; % element EI values w_e = [-2400. 0. 0.] ; % element line load % zero EBC and free DOF vector subscripts Fixed = [1 5 13 14] ; Free = [2 3 4 6 7 8 9 10 11 12]; n_c = max (size (Fixed)) ; % number of constraints % allocate required arrays S = zeros (n_d, n_d); c = zeros (n_d, 1); % system arrays S_e = zeros (n_i, n_i); c_f = zeros (n_i, 1); % elem arrays index = zeros (n_i, 1); % element scatter numbers u = zeros (n_d, 1); % system deflections and slopes React = zeros (n_c, 1); % system reactions % insert single point load and scatter it P_n = 5 ; % node of vertical point load P_eq = n_g*(P_n - 1) + 1 ; % eq of vertical load/displacement P_v = -1e4 ; % vertical point load value c(P_eq) = P_v ; % insert (scatter) point load into system for k = 1:n_e ; % begin element loop e=1 ====> ====> e_nodes = nodes (k, 1:3) ; % nodes on the element L = abs (x(e_nodes(3)) - x(e_nodes(1))) ; % compute length EI = EIs (k) ; % constant properties Line_Load = w_e(k) * [1 1 1]' ; % constant line load [index] = get_element_index(n_g, n_n, e_nodes);% sys numbers [S_e] = matrix_S_E_L3_C1 (EI, L) ; % bending stiffness [c_f] = matrix_c_f_L3_C1 (L, Line_Load) ; % load resultant S (index, index) = S (index, index) + S_e ; % S assembled c (index) = c (index) + c_f ; % c assembled end ; % for all k elements <==== <==== <==== <==== <==== fprintf ('\n') fprintf('Total transverse forces = %10.4e \n',sum(c(1:2:n_d))) % Save system rows for reaction recovery after solution RS = S(Fixed, :) ; RC = c(Fixed) ; % optional save method % independent equations K * U = C. Partition and solve K = S (Free, Free) ; C = c (Free) ; % known arrays U = K \ C ; % non-zero displacements u (Free) = U ; % gather back into full list list_save_results (n_g, n_m, u) ; % pretty print % graph_L3_C1_result (u, x, nodes) ; % plot the solution graph_L3_C1_moment (u, x, nodes, EIs) ; % plot the moment graph_L3_C1_shear (u, x, nodes, EIs) ; % plot the shear React = RS * u - RC ; % Reaction forces or moments Totals = zeros (n_g, 1) ; % zero system reaction components fprintf ('\n'); fprintf ('System Reactions \n') ; % header fprintf ('Node DOF Value \n') ; % header % Output reaction node_number, component_number, value for j = 1:n_c ; % extract all EBC reactions eq = Fixed(j) ; % equation number at reaction node = ceil(eq/n_g) ; % node at DOF j j_g = eq - (node - 1)*n_g ; % 1 <= j_g <= n_g value = React (j) ; % reaction value fprintf('%i %i %10.4e \n',node, j_g, value); % print Totals (j_g) = Totals (j_g) + value ; % sum all components end ; % for over all reaction values fprintf ('\n') fprintf('Total reaction transverse forces: %10.4e \n', ... sum(Totals(1))) fprintf('Total reaction couples: %10.4e \n', ... sum(Totals(2))) fprintf ('\n') % post process for element reactions fprintf('Element Reaction Values (none at mid-node) \n')% head fprintf('Elem Node Force Moment \n') % head for k = 1:n_e ; % begin element loop e=1 ====> ====> e_nodes = nodes (k, 1:3) ; % nodes on the element [index] = get_element_index(n_g, n_n, e_nodes);% sys numbers u_e = u(index) ; % answers for this element L = abs(x(e_nodes(3))-x(e_nodes(1))); % compute length % Recover or re-compute (here) the element matrices EI = EIs (k) ; % constant properties Line_Load = w_e(k) * [1 1 1]' ; % constant line load [S_e] = matrix_S_E_L3_C1 (EI, L) ; % bending stiffness [c_f] = matrix_c_f_L3_C1 (L, Line_Load) ; % load resultant % Assign any point source to the first element with that node for eq = 1:n_i ; % loop over element nodes n = index (eq) ; % find system DOF number if ( n == P_eq ) ; % found the input point value c_f (eq) = c_f (eq) + P_v ; % copy to this node end ; % if first appearance of this point source end ; % for node eq of the current element c_r = S_e * u_e - c_f ; % element reaction values format short e % make pretty print fprintf('%i %i %11.4e %11.4e \n', ... k, e_nodes(1), c_r(1), c_r(2)) % el node reactions fprintf('%i %i %11.4e %11.4e \n', ... k, e_nodes(3), c_r(5), c_r(6)) % el node reactions end ; % for all k elements <=== <==== <==== <==== <==== <==== % end Two_span_three_L3_C1_beam_react ======================== % Running gives: % Total transverse forces = -3.4000e+04 % % Computed Solution: % Node, 2 results per node % 1, 0.000e+00 -6.596e+04 % 2, -1.848e+05 8.511e+03 % 3, 0.000e+00 3.191e+04 % 4, -6.156e+04 -3.333e+04 % 5, -1.708e+05 1.603e+04 % 6, -7.337e+04 3.869e+04 % 7, 0.000e+00 0.000e+00 % % System Reactions % Node DOF Value % 1 1 9.9574e+03 % 3 1 1.7226e+04 % 7 1 6.8164e+03 % 7 2 -2.3121e+04 % % w=-2,400 lb/ft P=-10,000 lb % (1) (2) | (3) % wwwwwwwwwwwww V % *-----*-----*--------*--------*---*---* % 1 2 3 4 5 6 7 % ^ ^ ^ % | 9,957 | 17,226 lb 6,816 | / % ft-lb -23,121 <\-/ % % Total reaction transverse forces: 3.4000e+04 % Total reaction couples: -2.3121e+04 % % Element Reaction Values (none at mid-node) % Elem Node Force Moment % 1 1 9.9574e+03 -1.3642e-12 % 1 3 1.4043e+04 -2.0426e+04 % % w=2,400 lb/ft % (1) (24,000 lb total) % wwwwwwwwwwwww % *-----*-----* % ^ ^ ^ % | 9,957 | 14,043 lb % 0 -20426 ^\_/ ft-lb (counter clockwise) % % Elem Node Force Moment % 2 3 3.1836e+03 2.0426e+04 % 2 5 6.8164e+03 1.7778e+04 % % P=10,000 lb % (2) | % (17,226-14,043) V % *--------*--------* % ^ ^ % | 3,183.6 lb | 6,816.4 % +20426 \_/^ ft-lb +17,778 \-/^ % % Elem Node Force Moment % 3 5 3.1836e+03 -1.7778e+04 % 3 7 6.8164e+03 -2.3121e+04 % % (3) % (-10,000+6,816.4) % *---*---* % ^ ^ % 3,183.6 | | 6,816.4 % <\-/ -17,778 <\-/ -23,121 function [S_e] = matrix_S_E_L3_C1 (EI, L) %=================== % E = elastic modulus, I = moment of inertial, L = elem length L2 = L^2 ; L3 = L^3 ; % shorthand % Constant property quintic element bending stiffness S_e =[5092, 1138*L, -3584, 1920*L, -1508, 242*L ; 1138*L, 332*L2, -896*L, 320*L2, -242*L, 38*L2 ; -3584, -896*L, 7168, 0, -3584, 896*L ; 1920*L, 320*L2, 0, 1280*L2, -1920*L, 320*L2 ; -1508, -242*L, -3584, -1920*L, 5092, -1138*L ; 242*L, 38*L2, 896*L, 320*L2, -1138*L, 332*L2 ] ... *EI/(35 * L^3) ; % stiffness % end matrix_S_E_L3_C1 %==================================== function [c_f] = matrix_c_f_L3_C1 (L, Line_Load) % Quintic beam resultant shears and moments from % a quadratic line load through [f1 f2 f3] % L = element length % c_f = resultant at nodes, R_e * Line_Load % R_e = rectangular comversion matrix % Line_Load = [f1 f2 f3]' the load per length at % left, middle, and right node R_e = L * [57 44 -3 ; % 6 by 3 rectangle 3*L 4*L 0 ; 16 192 16 ; -8*L 0 8*L ; -3 44 57 ; 0 -4*L -3*L ] / 420 ; c_f = R_e * Line_Load ; % 6 by 1 resultant % end matrix_c_f_L3_C1 %========================