% function [] = Three_Layer_Wall () % revised 1/28/19 % Chimney wall with three brick layersmodelled with % three 2-node elements with sequential nodes % J.E. Akin, Rice University, Mech 417 % (K*A u,x),x + G*A = 0 for constant property elements % u1 u2 u3 u4 temperatures: % *-----------*-----------*-----------* ---> x x=x1 + L * r % 1 2 3 4 nodes % (1) (2) (3) elements % <---L1=9"---><--L2=5"--><---L3=7.5"-> lengths % K1=0.72 K2=0.08 K3=0.5 BTU/(hr ft F)) % A=1 ft^2 % 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 A = 1 ; % (ft^2) area G = 0. ; % BTU/(hr ft^3) gen rate per vol L = [9 5 7.5] / 12 ; % lengths , (ft) K = [0.72 0.08 0.5] ; % conductivity BTU/(hr ft F) c_P = [0 0 0 0]' ; % point heat flow (BTU/hr) % set given EBC locations and free locations Fix = [1 4] ; % all known EBC Free = [2 3] ; % active dof u(Fix) = [1500. 150.] ; % set temperature values for El = 1:n_e ; % Loop over all the elements ---> ---> ---> k_e = K(El)*A / L(El) ; % wall thermal stiffness S_e = k_e * [ 1, -1 ; % stiffness (BTU/hr F) -1, 1 ] ; % stiffness matrix c_e = G*A * L(El) * [1 1]' /2 ; % resultant source (BTU/hr) % 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 heat flows c' is \n");disp(c') fprintf ('Total applied source is %g (BTU/hr) \n', sum (c)) ; % Bring given temperature 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 temperature 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 temperatures u (Free) = u_f ; % update temperatures, for reactions fprintf ('\n') fprintf ('Node, Computed temperatures (F) \n') ; for k = 1:n_m ; % loop over all nodes --> --> fprintf ('%i, %10.4e \n', k, u(k)) ; end ; % for k nodes <-- <-- <-- x = [0 L(1) (L(1)+L(2)) (L(1)+L(2)+L(3))] graph_1d_result (n_e, n_m, n_n, x, nodes, u) % Recover the system reactions (at Fix nodes where EBC == 0) fprintf ('\n') ; % new line fprintf ('Reaction (+ entering, - exiting) heat flows (BTU) \n') n_c = max (size (Fix)) ; % number of EBC if ( n_c > 0 ) ; % thermal EBC given 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, heat flow %g (BTU) \n', Req, R (Req)) ; end ; % for k reactions <-- <-- <-- <-- fprintf ('Note: Total of all reactions is %g \n', sum (R)) end ; % if temperature given (vs. convection NBC) % Post-process for element temperature gradient and heat flux fprintf ('\n') ; % new line fprintf ('Element heat flux per area (BTU/ft^2) \n') for k = 1:n_e ; % Loop over all elements ===> ===> ===> ===> % Gather the temperatures for the element eqs = nodes (k, :) ; % system equation numbers u_e = u (eqs) ; % local temperatures % Temperature gradient, 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 Grad_T_e = B_e * u_e ; % constant gradient % Apply Fourier's law, s_e = -K_e * Grad_T_e s_e = -K(k) * Grad_T_e ; % heat flux per unit area fprintf ('Element %i heat flux is %10.4e \n', ... k, s_e) ; end ; % loop over elements <=== <=== <=== <=== <=== <=== % end Three_Layer_Wall ======================================= % Running gives % >> Three_Layer_Wall % After element 1 S and c' are % 0.9600 -0.9600 0 0 % -0.9600 0.9600 0 0 % 0 0 0 0 % 0 0 0 0 % % 0 0 0 0 % % After element 2 S and c' are % 0.9600 -0.9600 0 0 % -0.9600 1.1520 -0.1920 0 % 0 -0.1920 0.1920 0 % 0 0 0 0 % % 0 0 0 0 % % After element 3 S and c' are % 0.9600 -0.9600 0 0 % -0.9600 1.1520 -0.1920 0 % 0 -0.1920 0.9920 -0.8000 % 0 0 -0.8000 0.8000 % % 0 0 0 0 % % After adding external heat flows c' is % 0 0 0 0 % % Total applied source is 0 (BTU/hr) % After temperature BCs reduced S and c' are % 1.1520 -0.1920 % -0.1920 0.9920 % % 1440 120 % % % Node, Computed temperatures (F) % 1, 1.5000e+03 % 2, 1.3125e+03 % 3, 3.7500e+02 % 4, 1.5000e+02 % % Reaction (+ entering, - exiting) heat flows (BTU) % Node 1, heat flow 180 (BTU) % Node 4, heat flow -180 (BTU) % Note: Total of all reactions is 1.7053e-13 % % Element heat flux per area (BTU/ft^2) % Element 1 heat flux is 1.8000e+02 % Element 2 heat flux is 1.8000e+02 % Element 3 heat flux is 1.8000e+02