Why am I getting a error 'Dimensions of arrays being concatenated are not consistent'

2 views (last 30 days)
Why am i getting an error message. I thought i was producing a 12x12 zeros matrix, then adding the matrix calculcated to the zeroes. This works when I use a 4x4. why doesnt 12x12 work. LINE 117 is the error.
%
% BarFE.m
%
% The purpose of BarFE.m is to introduce you to coding your own
% Finite element software. Herein, we will use bar elements in 2D for
% simplicity.
%
% Use this file as a template to expand it to 3D, then further add
% DoF like shaft and beams.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Written by: Terence Macquart.
% Last Updated: 07/04/2021.
% Updated by:
% Change log:
%
%
% Contact details:
% BCI (ACCIS) - Bristol Composite Institute
% Department of Aerospace Engineering - University of Bristol
% Queen's Building - University Walk
% Bristol - BS8 1TR
% U.K.
% E-mail: terence.macquart@bristol.ac.uk
% Phone:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% swap folder and add any subfolders/files to your matlab path
if(~isdeployed), cd(fileparts(which('BarFE_Example')));end
Error using cd
Unable to change current folder to '' (Name is nonexistent or not a folder).
clear; close all; clc
maindir = pwd;
maindir(maindir=='\')='/'; % for linux
addpath (genpath(strcat(maindir)))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Step 1/2 geometry
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% node and element connectivity
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% node coordinates (set 3rd column to 0 for 2D)
nodes = [...
0, 0, 0 ;
1, 1, 0;
0, 1, 0;
1, 0, 0;
2, 1, 0;
2, 0, 0;
3, 1, 0];
% element connectivity matrix
eConn = [...
1, 2;
3, 2;
1, 4;
4, 5;
4, 6;
2, 5;
6, 7;
5, 7;
4, 2;
6, 5];
ne = size(eConn,1); % number of bar elements
% plot the structure nodes
figure('name','Nodes and Elements')
hold on
plot(nodes(:,1), nodes(:,2), 'blue o', 'markerfacecolor', 'blue')
xlabel('x (m)')
ylabel('y (m)')
axis equal
grid on
for i=1:size(nodes,1)
text(nodes(i,1)-0.05,nodes(i,2)-0.05,num2str(i),'color','red')
end
% compute theta and length for each element
Le = zeros(ne,1); % element length
theta = zeros(ne,1); % element orientation(rad)
for i=1:ne
elemtNode1 = nodes(eConn(i,1),:)';
elemtNode2 = nodes(eConn(i,2),:)';
opp = elemtNode2(2) - elemtNode1(2);
adj = elemtNode2(1) - elemtNode1(1);
theta(i) = atan2(opp,adj);
Le(i) = sqrt( sum((elemtNode2-elemtNode1).^2) ); % = norm(elemtNode2-elemtNode1)
% plot the structure elements
plot([elemtNode1(1), elemtNode2(1)], [elemtNode1(2), elemtNode2(2)],'black')
text(mean([elemtNode1(1), elemtNode2(1)]),mean([elemtNode1(2), elemtNode2(2)]),num2str(i))
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Step 3/4/5 Element stiffness matrices (bars only in this example)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% element properties and stiffness matrix
%%%%%%% all element are assumed to have the same area and Young's modulus
%%%%%%% Ke is the 2D element stiffness matrix in the element ref. frame
%%%%%%% Kg is the 2D element stiffness matrix in the global ref. frame
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A = 0.01; % bar cross-sectional area (m^2)
E = 190*1e9; % Young's modulus (Pa)
Izz = 5000;
Iyy = 5000;
G = 5000;
J = 5000;
Ke = zeros(12,12,ne);
Kg = zeros(12,12,ne);
for i=1:ne
Ke(:,:,i) = [A*E/(Le(i)), 0, 0, 0, 0, 0, -A*E/(Le(i)), 0, 0, 0, 0, 0;
0, 12*E*Izz/Le(i)^3, 0, 0, 0, 6*Le(i)*E*Izz/Le(i)^3, 0, -12*E*Izz/Le(i)^3, 0, 0, 0, 6*Le(i)*E*Izz/Le(i)^3;
0, 0, 12*E*Iyy/Le(i)^3, 0, -6*Le(i)*E*Iyy/Le(i)^3, 0, 0, 0, -12*E*Iyy/Le(i)^3, 0, -6*Le(i)*E*Iyy/Le(i)^3, 0;
0, 0, 0, G*J/Le(i), 0, 0, 0, 0, 0, -G*J/Le(i), 0, 0;
0, 0, -6*Le(i)*E*Iyy/Le(i)^3, 0, 4*Le(i)^2*E*Iyy/Le(i)^3, 0, 0, 0, 6*Le(i)*E*Iyy/Le(i)^3, 0, 2*Le(i)^2*E*Iyy/Le(i)^3, 0;
6*Le(i)*E*Izz/Le(i)^3, 0, 0, 0, 4*Le(i)^2*E*Izz/Le(i)^3, 0, -6*Le(i)*E*Izz/Le(i)^3, 0, 0, 0, 2*Le(i)^2*E*Izz/Le(i)^3;
A*E/Le(i), 0, 0, 0, 0, 0, A*E/Le(i), 0, 0, 0, 0, 0;
0, -12*E*Izz/Le(i)^3, 0 ,0, 0, -6*Le(i)*E*Izz/Le(i)^3, 0, 12*E*Izz/Le(i)^3, 0, -6*Le(i)*E*Izz/Le(i)^3;
0, 0, -12*E*Iyy/Le(i)^3, 0, 6*Le(i)*E*Iyy/Le(i)^3, 0, 0, 0, 12*E*Iyy/Le(i)^3, 0, 6*Le(i)*E*Iyy/Le(i)^3, 0;
0, 0, 0, -G*J/Le(i), 0, 0, 0, 0, 0, G*J/Le(i), 0, 0;
0, 0, -6*Le(i)*E*Iyy/Le(i)^3, 0, 2*Le(i)^2*E*Iyy/Le(i)^3, 0, 0, 0, 6*Le(i)*E*Iyy/Le(i)^3, 0, 4*Le(i)^2*E*Iyy/Le(i)^3, 0;
0, 6*Le(i)*E*Izz/Le(i)^3, 0, 0, 0, 2*Le(i)^2*E*Izz/Le(i)^3, 0, -6*Le(i)*E*Izz/Le(i)^3, 0, 0, 0, 4*Le(i)^2*E*Izz/Le(i)^3];
R = [cos(theta(i)) -sin(theta(i));
sin(theta(i)) cos(theta(i))];
R = blkdiag(R,R);
Kg(:,:,i) = R*Ke(:,:,i)*R';
end
%% Step 6 Assembly
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% element assembly
%%%%%%% the global stiffness matrix is 14x14 since we have 2 degrees of
%%%%%%% freedom per node (translations in x and y), and 7 nodes
%%%%%%%
%%%%%%% here I chose to number the degrees of freedom as follows
%%%%%%% Dof = description
%%%%%%% 1 = node 1 translation in x
%%%%%%% 2 = node 1 translation in y
%%%%%%% 3 = node 2 translation in x
%%%%%%% 4 = node 3 translation in y
%%%%%%% ...
%%%%%%% 14 = node 7 translation in y
%%%%%%%
%%%%%%% this choice is arbitrary and you could chose another one.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
K = zeros(14,14); % global stiffness matrix
for i=1:ne
nodeID = eConn(i,:); % index of nodes linked to element i
id = sort([nodeID*2-1 nodeID*2]); % index of dofs linked to element i
K(id,id) = K(id,id) + Kg(:,:,i); % add stiffnes of element i to global stiffness matrix
end
%% Step 7/8 Boundary conditions & Solve
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% Apply Boundary conditions & loads, then solve
%%%%%%% Clamped at node 1 and 3
%%%%%%% vertical downward force at node 7 Fy = -1e4 Newtons
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u = zeros(14,1);
F = zeros(14,1);
F(end) = -4e6; % external force (N) applied at the last dof
% remove the boundary conditions and solve the remaing system of equations
idBC = [1 2 5 6]; % dofs set to zero due to clamps
idFree = [3 4 7:14]; % remaining free dofs
Kf = K(idFree,idFree);
Ff = F(idFree);
u(idFree) = Kf^-1*Ff; % displacement solution caused by F
%% Step 9 Post-processing
% format the output and visualise the results
u_xy = reshape(u,2,[])';
DeformedNode = nodes(:,1:2)+u_xy;
L = zeros(ne,1); % deformed element length
for i=1:ne
elemtNode1 = DeformedNode(eConn(i,1),:)';
elemtNode2 = DeformedNode(eConn(i,2),:)';
L(i) = sqrt( sum((elemtNode2-elemtNode1).^2) );
if L(i)>Le(i)
% member in tension
plot([elemtNode1(1), elemtNode2(1)], [elemtNode1(2), elemtNode2(2)],'blue','linewidth',2)
else
% member in compression
plot([elemtNode1(1), elemtNode2(1)], [elemtNode1(2), elemtNode2(2)],'red','linewidth',2)
end
end
% strains due to axial deformation (from element summary table for bar)
eps_xx = zeros(ne,1);
for i=1:ne
% for each element find the nodes and DoFs indices
idNodes = eConn(i,:)';
Dof_Node1 = [idNodes(1)*2-1, idNodes(1)*2];
Dof_Node2 = [idNodes(2)*2-1, idNodes(2)*2];
ug = u([Dof_Node1,Dof_Node2]); % global displacement of the nodes
% convert ug into displacement of the nodes expressed in the element coordinate system
R = [cos(theta(i)) -sin(theta(i));
sin(theta(i)) cos(theta(i))];
R = blkdiag(R,R);
ue = R'*ug;
% evaluate the element axial strain
% you can check that eps_xx = (L(i)-Le(i))/Le(i) for small displacements
eps_xx(i) = [-1/Le(i) 1/Le(i)]*ue([1 3]);
end
%% Step 10 Validate
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%% Validation against Abaqus results
% Node Label U.U1 U.U2
% @Loc 1 @Loc 1
% -------------------------------------------------
% 1 -6.31579E-03 -37.1723E-03
% 2 -4.21053E-03 -14.3756E-03
% 3 6.31579E-03 -12.2704E-03
% 4 10.5263E-03 -35.0671E-03
% 5 -12.0000E-30 -4.00000E-30
% 6 12.6316E-03 -62.0743E-03
% 7 12.0000E-30 0.
%
% abaqus node numbering is different than our so we need to make sure to
% compare the correct nodes together --> this is why we use idmap
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
AbaqusData = [ 1 -6.31579E-03 -37.1723E-03
2 -4.21053E-03 -14.3756E-03
3 6.31579E-03 -12.2704E-03
4 10.5263E-03 -35.0671E-03
5 -12.0000E-30 -4.00000E-30
6 12.6316E-03 -62.0743E-03
7 12.0000E-30 0.];
idmap = [5 3 7 2 4 1 6];
error = abs(u_xy-AbaqusData(idmap,2:3));
assert(all(error(:)<1e-6),'benchmark AbaqusData not maching the results')

Answers (1)

Voss
Voss on 20 Nov 2022
This line has only 11 elements:
6*Le(i)*E*Izz/Le(i)^3, 0, 0, 0, 4*Le(i)^2*E*Izz/Le(i)^3, 0, -6*Le(i)*E*Izz/Le(i)^3, 0, 0, 0, 2*Le(i)^2*E*Izz/Le(i)^3;
And this line has only 10 elements:
0, -12*E*Izz/Le(i)^3, 0 ,0, 0, -6*Le(i)*E*Izz/Le(i)^3, 0, 12*E*Izz/Le(i)^3, 0, -6*Le(i)*E*Izz/Le(i)^3;
The other 10 lines each have 12 elements.

Categories

Find more on Polymers in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by