% function [] = Bar_pt_line_temp () % revised 2/24/2020 % [CORRECTED VERSION POSTED AS HW7 SOLUTION 3/20/20] % The 1D solution of [d/dx (E*A du/dx)] - w =0 % Bar with axial point, line, and or thermal loads % For type 1 stress element with nodes 1---2---3 -->r % (Element interfaces must occur at any discontinuity) % Input coordinates, any point loads, connection list % [Free hanging bar with line load w = rho*area] x = [0.0, 12., 24.] ; % all coords, (m) Pt = [0.0, 0.0, 0.0] ; % all external pt forces (N) nodes = [1, 2, 3] ; % mesh connections % Give the essential (Dirichlet) BC nodes and values EBC_n = [1, 3] ; % node at EBC EBC_v = [0, 0] ; % EBC value (m) % Set controls integers n_c = size (EBC_n, 2) ; % number of contraints n_C = size (EBC_v, 2) ; % number of contraints n_e = size (nodes, 1) ; % number of elements n_m = max ( max (nodes)) ; % number of nodes n_P = size (Pt, 2) ; % number of point loads n_n = size (nodes, 2) ; % number of nodes per el n_g = 1 ; % number of DOF per node n_p = 1 ; % dimension of parametric space n_r = 1 ; % number of rows in element B_e matrix n_s = 1 ; % dimensions of physical space if ( n_m ~= n_P ) ; % then user data error error ('Number of nodes and pt. forces do not match') end ; % if bad BC data if ( n_c ~= n_C ) ; % then user data error error ('Number of BC nodes and values do not match') end ; % if bad BC data % Compute derived controls n_d = n_g * n_m ; % system degrees of freedom (DOF) n_i = n_g * n_n ; % number of DOF per element % Set logics post = 1 ; % turn on(1)/off(0) post-processing option % Give all element's area, modulus, line_load, alpha, del_T % (m^2), (N/m^2), (N/m), (1/C), (C) e_A = [0.25] ; % (in^2) e_E = [29.e6] ; % (psi) elastic modulus e_w = [0.] ; % (lb/in) line load e_al = [6.60e-6] ; % (1/F) coeff. thermal expansion e_dT = [120-60] ; % (F) temp. change from stress free % Allocate (main) system memory Ans = zeros (n_d, 1) ; % system results C = zeros (n_d, 1) ; % system force or source S = zeros (n_d, n_d) ; % system stiffness R_C = zeros (n_c, 1) ; % reaction C partition R_S = zeros (n_c, n_d) ; % reaction S partition React = zeros (n_c, 1) ; % reaction values % Insert any point loads into system loads C (:) = Pt (:) ; % assemble pt loads Total_P = sum(Pt) ; % get load sum Total_w = 0 ; % zero load sum Total_T = 0 ; % zero load sum % Allocate and clear element type arrays c_e = zeros (n_i, 1) ; % clear array el sources S_e = zeros (n_i, n_i) ; % clear array el stiffness x_e = zeros (n_n, n_s) ; % element coordinates for k = 1:n_e ; % loop over elements ====>> ====>> ====>> % Gather coordinates, properties, form element arrays E_e = e_E(k) ; % element's modulus A_e = e_A(k) ; % element's area w_e = e_w(k) ; % element's line load al_e = e_al(k) ; % C.O.T.E. dT_e = e_dT(k) ; % change in temp. f_T = E_e*A_e*al_e*dT_e ; % thermal force rows = nodes (k, 1:n_n) ; % element's connections x_e = x (rows) ; % gather x-coords L_e = max (x_e) - min (x_e) ; % element length S_e = [ 7, -8, 1, ; ... % stiffness L3_C0 (N/m) -8, 16, -8, ; ... 1, -8, 7 ] * E_e * A_e / (3*L_e) ; % stiffness c_w = w_e * L_e * [1; 4; 1] / 6 ; % line load force (N) c_t = f_T * [-1; 0; 1] ; % thermal load (N) c_e = c_w + c_t ; % total el load (N) Total_w = Total_w + sum(c_w) ; % update line loads Total_T = Total_T + f_T ; % update thermal loads % Scatter element arrays into system arrays % rows = vector subscript converts el to system eq numbers S (rows, rows) = S (rows, rows) + S_e (:, :) ; % add stiff C (rows) = C (rows) + c_e (:) ; % add forces end ; % for each k element in mesh <<==== <<==== <<==== <<==== fprintf ('Total external point loads = %8.3e \n', Total_P) fprintf ('Total external line loads = %8.3e \n', Total_w) fprintf ('Total thermal loadings = %8.3e \n', Total_T) % ---- Enforce (Dirichlet) EBC here ---_ % by trick to avoid partitions, save reaction row first for k = 1: n_c ; % Loop over Dirichlet BC --> --> --> --> --> eq = EBC_n(k) ; val = EBC_v(k) ; % EBC equation and value R_S(k,:) = S (eq, 1:n_d) ; R_C(k) = C(eq) ; % for later use Diag = max ( diag(S) ) ; % for better conditioning % Carry known columns*EBC to RHS. Zero that column and row C (1:n_d) = C (1:n_d) - val * S (1:n_d, eq) ; % new RHS S (1:n_d, eq) = 0 ; S (eq, 1:n_d) = 0 ; % new LHS % Insert EBC identity, Diag * EBC_dof = Diag * val C (eq) = Diag * val ; % new RHS S (eq, eq) = Diag ; % new LHS end ; % for k linear constraint equations <-- <-- <-- <-- <-- % ----------------- Dirichlet BC enforced -------------------- % Solve the modified system Ans = S \ C ; % Results fprintf ('Node, x, solution \n') % title for k = 1:n_m % for all nodes % pretty x_p = x(k) ; % check fprintf ('%i, %7.2e, %8.3e \n', k, x(k), Ans(k)) % print end ; % for all k system nodes (and 1 DOF) % Get reaction (NBC) at EBC node (via trick method) React = R_S * Ans - R_C ; % NBC value fprintf ('\n') for k = 1:n_c ; % list reactions fprintf('Reaction at node %i = %8.3e \n',EBC_n(k),React(k, 1)) end ; % for k restraints % post-process the results ?? if ( post == 1 ) % use computed solution p_n = 2 ; % number of post-process points in element r_q (1) = 2.1132486540518711774543e-01 ; % from tables r_q (2) = 7.8867513459481288225457e-01 ; % from tables fprintf ('\n') fprintf ('Post-processing Results (at 2 qp) \n') fprintf ('Elem x strain stress force \n') % Allocate and clear element type arrays x_e = zeros (n_n, n_s) ; % element coordinates Ans_e = zeros (n_i, 1) ; % element answers for k = 1:n_e ; % loop over elements ====>> ====>> ====>> % Gather coordinates, properties, form element arrays rows = nodes (k, 1:n_n) ; % this element connections x_e = x (rows) ; % gather x-coords L_e = max (x_e) - min (x_e) ; % element length Ans_e = Ans (rows) ; % gather element solution % Loop over most accurate points (quadrature) & node points for q = 1:p_n ; % loop over points ---> ---> ---> ---> ---> r = r_q(q) ; % non-dimensional pt % Find interpolations, local derivatives, and location H_q = [(1-3*r+2*r^2),(4*r-4*r^2),(2*r^2-r)] ; % L3 C0 DLH_q = [(-3+4*r), (4-8*r), (4*r-1)] ; % L3 C0 x_q = H_q * x_e' ; % physical x (m) % Find Jacobian, inverse, physical derivative Jac = DLH_q * x_e' ; % Jacobian inv_J = 1 / Jac ; % inverse Jacobian dudx_q = inv_J * DLH_q * Ans_e ; % gradient (strain m/m) % ******************************************************* strain_m = dudx_q ; % mechanical strain strain_T = e_al(k) * e_dT(k) ; % thermal strain strain_n = strain_m - strain_T ; % net strain % ******************************************************* stress_q = E_e * strain_n ; % stress (N/m^2) force_q = A_e * stress_q ; % force (N) fprintf ('%i, %7.2e %8.3e %8.3e %8.3e \n', ... k, x_q, strain_n, stress_q, force_q) % FEA values end ; % for q selected points <--- <--- <--- <--- <--- <--- end ; % for each k element in mesh <<==== <<==== <<==== <<==== fprintf ('\n') fprintf ('The above force answers are sometimes NOT consistent \n') fprintf ('with the correct reactions above due to a program bug \n') fprintf ('found 3/6/20 and to be announced. \n') fprintf ('Can you find the bug(s) first ?? \n') end ; % if use computed solution % end Bar_pt_line_temp % Running gives %=================================================== % Hanging bar, top fixed, no point or thermal loads % x = [0.0, 2.0, 4.0] ; % coords, (m) % Pt = [0.0, 0.0, 0.0] ; % external force (N) % nodes = [1, 2, 3] ; % mesh connections % EBC_n = [1] ; % node at EBC % EBC_v = [0] ; % EBC value (m) % e_A = [2.0] % e_E = [10.] % e_w = [3.] % e_al = [0.] % e_dT = [0.] % % Total external point loads = 0.000e+00 % Total external line loads = 1.200e+01 % Total thermal loadings = 0.000e+00 % Node, x, solution % % 1, 0.00e+00, 0.000e+00 % 2, 2.00e+00, 9.000e-01 % 3, 4.00e+00, 1.200e+00 % % Reaction at node 1 = -1.200e+01 % % Post-processing Results (at 2 qp) % Elem x strain stress force % 1, 8.45e-01 4.732e-01 4.732e+00 9.464e+00 % 1, 3.15e+00 1.268e-01 1.268e+00 2.536e+00 % %=================================================== % Repeat with both ends fixed % x = [0.0, 2.0, 4.0] ; % coords, (m) % Pt = [0.0, 0.0, 0.0] ; % external force (N) % nodes = [1, 2, 3] ; % mesh connections % EBC_n = [1, 3] % EBC_v = [0, 0] ; % two EBC % e_A = [2.0] % e_E = [10.] % e_w = [3.] % e_al = [0.] % e_dT = [0.] % % Total external point loads = 0.000e+00 % Total external line loads = 1.200e+01 % Total thermal loadings = 0.000e+00 % % Node, x, solution % 1, 0.00e+00, 0.000e+00 % 2, 2.00e+00, 3.000e-01 % 3, 4.00e+00, 0.000e+00 % % Reaction at node 1 = -6.000e+00 % Reaction at node 3 = -6.000e+00 % ******************************************************* % Post-processing Results (at 2 qp) % Elem x strain stress force % 1, 8.45e-01 1.732e-01 1.732e+00 3.464e+00 % 1, 3.15e+00 -1.732e-01 -1.732e+00 -3.464e+00 % %=================================================== % Repeat both ends fixed, uniform temperature rise % Hibbler, 'Mech. of Materials', 4th, Ex 4-10 % x = [0.0, 12., 24.] ; % all coords, (m) % Pt = [0.0, 0.0, 0.0] ; % all external pt forces (N) % nodes = [1, 2, 3] ; % mesh connections % EBC_n = [1, 3] ; % node at EBC % EBC_v = [0, 0] ; % EBC value (m) % e_A = [0.25] ; % (in^2) area % e_E = [29.e6] ; % (psi) elastic modulus % e_w = [0.] ; % (lb/in) line load % e_al = [6.60e-6] ; % (1/F) coeff. thermal expansion % e_dT = [120-60] ; % (F) temp. change from stress free % % Total external point loads = 0.000e+00 % Total external line loads = 0.000e+00 % Total thermal loadings = 2.871e+03 % % Node, x, solution % 1, 0.00e+00, 0.000e+00 % 2, 1.20e+01, 0.000e+00 % 3, 2.40e+01, 0.000e+00 % % Reaction at node 1 = 2.871e+03 % Reaction at node 3 = -2.871e+03 % % Post-processing Results (at 2 qp) % Elem x strain stress force % 1, 5.07e+00 -3.960e-04 -1.148e+04 -2.871e+03 % 1, 1.89e+01 -3.960e-04 -1.148e+04 -2.871e+03 %===================================================