syms E A I L theta_t theta_b Fy6 real
%assumes all the beams are the same, can be changed.
beam_config = [1, 3, pi/2+theta_t
               2, 3, pi/2-theta_t
               3, 4, pi
               3, 5, 0
               4, 6, 0
               5, 6, pi
               6, 7, pi
               6, 8, 0
               7, 9, 0
               8, 9, pi
               9, 10, -pi/2-theta_b
               9, 11, -pi/2+theta_b];
K_actual = assemble_stiffness(11, beam_config, E, A, I, L);% Define stiffness matrix

constrained_nodes = [1, 2, 10, 11]; % Nodes that are fixed
constrained_dofs = [];
for node = constrained_nodes
    constrained_dofs = [constrained_dofs, 3*node-2, 3*node-1, 3*node]; % x, y, theta DOFs
end

% Define free DOFs (all DOFs minus constrained DOFs)
all_dofs = 1:size(K_actual, 1);
free_dofs = setdiff(all_dofs, constrained_dofs);

%Define applied Force
F = sym(zeros(33, 1)); % Only F17 is nonzero, y6 direction
F(17) = Fy6;

% Reduce stiffness matrix and force vector
K_reduced = K_actual(free_dofs, free_dofs) % Remove constrained rows/columns
F_reduced = F(free_dofs); % Remove constrained rows


% Solve the reduced system
%%%%%%%%%%%%%%%%%%%%%%%%%%%% GETS STUCK HERE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x_reduced = K_reduced \ F_reduced; % Solve for displacements at free DOFs


% Reconstruct the full displacement vector
x_full = zeros(size(K_actual, 1), 1);
x_full(free_dofs) = x_reduced; % Fill in free DOFs (constrained DOFs remain zero)

% Compute effective stiffness
K_eff = simplify(Fy6 / x_full(17)); % y6 corresponds to DOF 17
disp('Effective stiffness:');
disp(K_eff);

clear; clc;
syms ka kb kc kd ke F5
% 5 nodes, example parallel series config:
% /| |-->2-->|
% /|1|       |4-->5
% /| |-->3-->|
%
% node 1 is constrained.
% SOLVE FOR EFFECTIVE X STIFFNESS AT NODE 5
%
% test solution from series-parallel equations:
% K_test = inv(inv(inv(inv(ka)+inv(kb))+inv(inv(kc)+inv(kd)))+inv(ke))

% Define stiffness matrix
K = [ka+kc     -ka     -kc     0    0
     -ka    ka+kb       0   -kb     0
      -kc     0     kc+kd  -kd      0
      0       -kb   -kd    kb+kd+ke -ke
      0       0     0    -ke  ke]

% Define force vector
F = [0; 0; 0; 0; F5]; % Only F5 is nonzero

% Define constraint matrix (C)
% We enforce x1 = 0, so the constraint row is [1 0 0 0]
C = [1 0 0 0 0];
N_c = size(C,1)
% Define augmented matrix A and augmented RHS vector b
A = [K, C.'; C, zeros(N_c, N_c)]
b = [F; zeros(N_c,1)]; % 0 corresponds to the fixed value of x1

% Solve the augmented system
x = A \ b;

% Extract x5 (displacement at node 5)
x5 = x(5);

% Compute effective stiffness
K_eff = simplify(F5 / x5);

disp('Effective stiffness K_eff:');
disp(K_eff);
K_test = inv(inv(inv(inv(ka)+inv(kb))+inv(inv(kc)+inv(kd)))+inv(ke))%from series-parallel equations
isAlways(K_eff== K_test) % are they the same?

function K_global = assemble_stiffness(num_nodes, beam_data, E, A, I, L)
    % Inputs:
    % num_nodes - Total number of nodes in the structure
    % beam_data - Matrix where each row is [node1, node2, angle (radians)]
    % E - Young's modulus
    % A - Cross-sectional area
    % I - Second moment of area
    % L - Length of each beam (assumed the same for simplicity)

    % Initialize global stiffness matrix
    K_global = sym(zeros(3 * num_nodes));

    % Loop through each beam
    for i = 1:size(beam_data, 1)
        % Extract beam data
        node1 = beam_data(i, 1);
        node2 = beam_data(i, 2);
        angle = beam_data(i, 3);

        % Compute local stiffness matrix for the beam
        K_local = beam_stiffness_matrix(E, A, I, L);

        % Compute rotation matrix
        R = rotation_matrix(angle);

        % Transform local stiffness matrix to global coordinates
        K_global_beam = R' * K_local * R;

        % Assemble into the global stiffness matrix
        dof_indices = [3*node1-2, 3*node1-1, 3*node1, 3*node2-2, 3*node2-1, 3*node2];
        for ii = 1:6
            for jj = 1:6
                K_global(dof_indices(ii), dof_indices(jj)) = ...
                    K_global(dof_indices(ii), dof_indices(jj)) + K_global_beam(ii, jj);
            end
        end
    end
end

function K_local = beam_stiffness_matrix(E, A, I, L)
    % Returns the 6x6 local stiffness matrix for a beam element
    K_local = sym([
        E*A/L,        0,            0,         -E*A/L,        0,          0;
        0,     12*E*I/L^3,   6*E*I/L^2,        0,   -12*E*I/L^3,   6*E*I/L^2;
        0,      6*E*I/L^2,    4*E*I/L,         0,    -6*E*I/L^2,    2*E*I/L;
        -E*A/L,        0,            0,         E*A/L,         0,          0;
        0,    -12*E*I/L^3,  -6*E*I/L^2,        0,    12*E*I/L^3,  -6*E*I/L^2;
        0,      6*E*I/L^2,    2*E*I/L,         0,    -6*E*I/L^2,    4*E*I/L;
    ]);
end

function R = rotation_matrix(angle)
    % Returns the 6x6 rotation matrix for a beam rotated by `angle`
    c = cos(angle);
    s = sin(angle);

    % Rotation matrix for 2D with moments
    R = [
        c,  s,  0,  0,  0,  0;
       -s,  c,  0,  0,  0,  0;
        0,  0,  1,  0,  0,  0;
        0,  0,  0,  c,  s,  0;
        0,  0,  0, -s,  c,  0;
        0,  0,  0,  0,  0,  1;
    ];
end

Published with MATLAB® R2023b