function [] = DC_wires () % Ten wire DC Network by FEA % J.E. Akin, Rice University revised 1/29/2020 % unknowns are voltages, reactions are DC currents % Define element connectivity (like msh_typ_nodes.txt) % (1) (3) El Nodes R % 120 v BC *--k2--*--k6--* (6) 1 1-2 36. % | | | 2 1-3 1. % | k4 k9 3 2-5 2. % | | | 4 3-4 10. % k1 (4)*--k7--* (7) 5 4-5 20. % | | | 6 3-6 3. % | k5 k10 7 4-7 5. % | | | 8 5-8 4. % 0 v BC *--k3--*--k8--* (8) 9 6-7 7. % (2) (5) 10 7-8 6, % (j) represents wire connect % Define wire ends 'connection list' (line msh_typ_nodes.txt) nodes = [1 2; 1 3; 2 5; 3 4; 4 5; 3 6; 4 7; 5 8; 6 7; 7 8] % get the number of elements, nodes, and number of equations n_e = size (nodes, 1) ; % number of elements in mesh n_m = max (max(nodes)) ; % number of nodes in mesh n_g = 1 ; % number of unknowns per node n_d = n_g * n_m ; % nunmber of equations in the system fprintf ('System has %i elements \n', n_e) fprintf ('System has %i nodes \n', n_m) fprintf ('System has %i unknowns per node \n', n_g) fprintf ('System has %i equations \n', n_d) % Define element properties (like msh_properties.txt) % fprintf ('Wire resistances \n') R_j = [36. 1. 2. 10. 20. 3. 5. 4. 7. 6.] % Resistances % Define fixed and free voltage degrees of freedom (dof) % via equation numbers and given BC values (like msh_ebc.txt) Fixed = [1 2]; Free = [3 4 5 6 7 8] ; % (for n_g=1) BC_is = [120. 0.] ; % any boundary condition values at Fixed % *** % (Edit above two lines to change BC locations or values) % *** fprintf ('Given DC voltages at nodes \n'); disp(Fixed) fprintf ('are \n'); disp(BC_is) % Always allocate number of voltage, U, forces, F, and % the square matrix size K in K U = F assembly U = zeros (n_d, 1); F = zeros (n_d, 1); K = zeros (n_d, n_d); % Set the known voltage values (copy via vector subscript) U(Fixed) = BC_is (:) ; % (m) copied into BC locations % Set non-zero known external currents (like load_pt.txt) Applied = sum(F) ; % later compare to reactions fprintf ('Sum of external currents = %10.4e \n', Applied) % Loop over elements to assemble (singular) K matrix for n = 1:n_e ; % ====> ====> ====> ====> ====> ====> ====> % Gather element properties (and usually its coordinates) el_nodes = nodes (n, :) ; % gather element connections k = 1/R_j (n) ; % gather element properties ; % (Usually a numerical integration here.) K_e = k * [1 -1; % compute el stiffness -1 1] ; % compute el stiffness % Scatter element stiffness to system stiffness K (el_nodes, el_nodes) = K (el_nodes, el_nodes) + K_e ; end ; % for n-th element <==== <==== <==== <==== <==== <==== % These matrix equilibrium equations are not unique % because they do not yet include Dirichlet BCs.) fprintf('Assembled (singular) conductance matrix \n'); disp(K) fprintf('External input sources \n'); disp(F) % Transfer known BC voltage effects to Free loads F (Free) = F (Free) - K (Free, Fixed)*U (Fixed) % List the non-singular partition of the matrix system fprintf ('Free conductance partition \n'); disp(K(Free,Free)) fprintf ('Free currents \n'); disp(F(Free)) % Invert non-singular K portion, multiply it times the % Free forces to find the Free voltage U (Free) = K (Free, Free) \ F (Free) ; % (m) fprintf ('System DC voltages: \n') ; disp(U) % Optionally find system reactions at the Fixed voltage % needed to support the Dirichlet BCs F(Fixed) = K(Fixed, Fixed)*U(Fixed) + K(Fixed, Free)*U(Free) ; fprintf ('System currents at Fixed BC: \n'); disp(F(Fixed)) fprintf ('System currents sum = %10.4e \n', sum(F(Fixed))); fprintf ('Sum of external currents = %10.4e \n', Applied) %% Always find the element gradients, here Delta_U/L %% Loop over elements to compute the gradient %for n = 1:n_e ; % ====> ====> ====> ====> ====> ====> ====> % el_nodes = nodes (n, :) ; % gather element nodes % L_e = L_j (n) ; % gather element size % U_e = U(el_nodes) ; % gather end voltage % g_e = (U_e(2) - U_e(1)) / L_e ; % element gradient % fprintf ('Element %i gradient = %10.4e \n', n, g_e) %end ; % for n-th element <==== <==== <==== <==== <==== <==== % Optionally find the wire end currents (el reactions) % Loop over elements to compute their reactions (for n_g=1) fprintf ('Current entering (+) or leaving (-) end nodes \n') for n = 1:n_e ; % ====> ====> ====> ====> ====> ====> ====> el_nodes = nodes (n, :) ; % gather element connections U_e = U(el_nodes) ; % gather the end voltages k = 1/R_j (n) ; % gather element properties K_e = k * [1 -1; -1 1] ; % compute el stiffness R_e = K_e * U_e ; % recover end reaction currents fprintf ('Currents at element number %i \n', n) fprintf ('at nodes %i %i \n', nodes (n, :)) disp(R_e) end ; % for n-th element <==== <==== <==== <==== <==== <==== % Running gives % System currents at Fixed BC: % 11.4608 % -11.4608 % end DC_wires