% function [] = Axial_three_L2 () % revised 1/23/19 % bar modelled with three 2-node elements with sequential nodes % J.E. Akin, Rice University, Mech 417 % (E*A u,x),x + G*A = 0 for constant property elements % u1 u2 u3 u4 displacements % *-----------*-----------*-----------* ---> x x=x1 + L * r % 1 2 3 4 nodes % (1) (2) (3) elements % <---L1=2'--><---L2=3'--><---L3=4'---> lengths % % data to calculate storage and controls n_e = 3 ; % number of elements for this application n_g = 1 ; % dof per node (C0) n_n = 2 ; % number of element nodes (L2) n_m = 4 ; % number of mesh nodes n_i = n_g * n_n ; % dof per element n_d = n_g * n_m ; % dof in the full system % allocate element array sizes nodes = zeros (n_e, n_n) ; % element connections S_e = zeros (n_i, n_i) ; c_e = zeros (n_i, 1) ; u_e = zeros (n_i, 1) ; % element arrays % allocate system array sizes (before EBC) S = zeros (n_d, n_d) ; c = zeros (n_d, 1) ; % system arrays u = zeros (n_d, 1); R = zeros (n_d, 1); c_P = zeros (n_d, 1); % element connectivity data & constant properties nodes = [1 2 ; % element 1 connection list 2 3 ; % element 2 connection list 3 4 ] ; % element 3 connection list E = 30e6 ; % (psi) elastic modulus A = 1 ; % (in^2) area G = 0. ; % (lb/in^3) wt per volume` L = 12 * [2.0 3.0 4.0] ; % lengths , (inch) c_P = [-10000 3000 -2000 9000]'; % external point loads (lb) % set zero EBC locations and free locations Fix = [1] ; % all zero EBC Free = [2 3 4] ; % active dof u(Fix) = 2. ; % set displacement value for El = 1:n_e ; % Loop over all the elements ---> ---> ---> k_e = E*A / L(El) ; % bar axial stiffness S_e = k_e * [ 1, -1 ; % stiffness L2_C0 -1, 1 ] ; % stiffness matrix c_e = G*A * L(El) * [1 1]' /2 ; % resultant source % Assemble into equilibrium: S * u = c_P + c_G + c_BC e_nodes = nodes (El, :) ; % connected nodes Eqs = e_nodes ; % n_g = 1 only c (Eqs) = c (Eqs) + c_e ; % scatter sources S (Eqs, Eqs) = S (Eqs, Eqs) + S_e ; % scatter stiffness % Summary of element scatters fprintf ("After element %i S and c' are \n", El) disp(S) ; disp(c') end % for all elements <--- <--- <--- <--- <--- <--- <--- <--- c = c + c_P ; % add point loads fprintf ("After adding external point loads c' is \n");disp(c') fprintf ('Note: Total applied source is %g \n', sum (c)) ; % Bring given displacement BC to right hand side c (Free) = c (Free) - S (Free, Fix) * u (Fix) ; % Extract the independent matrix partitions S_f = S (Free, Free) ; % extract free independent stiffness c_f = c (Free) ; % extract independent forces % Reduced (non-singular) matrix equilibrium fprintf ("After displacement BCs reduced S and c' are \n") disp(S_f) ; disp(c_f') % Solve for free dof (omitting known value essential BC) u_f = S_f \ c_f ; % solve for independent displacements u (Free) = u_f ; % update displacements, for reactions fprintf ('\n') fprintf ('Node, Computed solution \n') ; %B format long for K = 1:n_m ; % loop over all nodes --> --> fprintf ('%i, %10.4e \n', K, u(K)) ; end ; % for K nodes <-- <-- <-- % Recover the system reactions (at Fix nodes where EBC == 0) fprintf ('\n') ; % new line n_c = max (size (Fix)) ; % number of EBC for k = 1:n_c ; % loop over number of EBC --> --> --> Req = Fix (k) ; % equation # R (Req) = S (Req, :) * u (:) - c (Req) ; % reaction value fprintf ('Node %i, Reaction %g \n', Req, R (Req)) ; end ; % for k reactions <-- <-- <-- <-- fprintf ('Note: Total of all reactions is %g \n', sum (R)) % Post-process for element strains and stresses fprintf ('\n') ; % new line for k = 1:n_e ; % Loop over all elements ===> ===> ===> ===> % Gather the displacements for the element eqs = nodes (k, :) ; % system equation numbers u_e = u (eqs) ; % local displacements % Mechanical strain, em = B_e * u_e, B_e = [-1 1]/L_e L_e = L (k) ; % element length B_e = [-1 1]/L_e ; % a constant em_e = B_e * u_e ; % constant strain % Apply stress-strain law, s_e = E_e * (em -et)_e s_e = E * (em_e - 0.0) ; % no thermal strain, et fprintf ('Element %i axial stress is %10.4e \n', ... k, s_e) ; end ; % loop over elements <=== <=== <=== <=== <=== <=== % end Axial_three_L2 ======================================= % Running gives % >> Axial_three_L2 % After element 1 S and c' are % 1250000 -1250000 0 0 % -1250000 1250000 0 0 % 0 0 0 0 % 0 0 0 0 % % 0 0 0 0 % % After element 2 S and c' are 1.0e+06 * % 1.2500 -1.2500 0 0 % -1.2500 2.0833 -0.8333 0 % 0 -0.8333 0.8333 0 % 0 0 0 0 % % 0 0 0 0 % % After element 3 S and c' are 1.0e+06 * % 1.2500 -1.2500 0 0 % -1.2500 2.0833 -0.8333 0 % 0 -0.8333 1.4583 -0.6250 % 0 0 -0.6250 0.6250 % % 0 0 0 0 % % After adding external point loads c' is % -10000 3000 -2000 9000 % % Note: Total applied source is 0 % % After BCs reduced S and c' are 1.0e+06 * % 2.0833 -0.8333 0 % -0.8333 1.4583 -0.6250 % 0 -0.6250 0.6250 % % 2503000 -2000 9000 % % % Node, Computed solution % 1, 2.0000e+00 % 2, 2.0080e+00 % 3, 2.0164e+00 % 4, 2.0308e+00 % % Node 1, Reaction -4.65661e-10 % Note: Total of all reactions is -4.65661e-10 % % Element 1 axial stress is 1.0000e+04 % Element 2 axial stress is 7.0000e+03 % Element 3 axial stress is 9.0000e+03