Index in position 1 exceeds array bounds (must not exceed 8).

2 次查看（过去 30 天）
Syed
Syed 2023-12-26
编辑： Dyuman Joshi 2023-12-26
Ran in:
I am trying to find the global stiffness matrix of 4 Q4 elements, not sure why im getting this error can someone help, kinda very noob at coding
this is my code:
%(a)
% Material properties
E = 207e9; % Young's modulus in Pa
nu = 0.32; % Poisson's ratio
% Pressure applied on the plate
P = 39 * 1e6; % Pa
% Element properties
t = 20*10^-3; % Plate thickness in meters
% Gauss quadrature points and weights
dof = 2; % Number of integration points in each direction
xiVec = [-1/sqrt(3), 1/sqrt(3)]; % Gauss quadrature points in the xi direction
etaVec = xiVec; % Gauss quadrature points in the eta direction
wVec = [1, 1]; % Gauss weights in both directions
% Initialize the element stiffness matrices
Ke1 = zeros(8, 8);
Ke2 = zeros(8, 8);
Ke3 = zeros(8, 8);
Ke4 = zeros(8, 8);
% Loop over the integration points
for i = 1:dof
for j = 1:dof
xi = xiVec(i);
eta = etaVec(j);
% Shape functions and derivatives
N = 1/4 * [(1 - xi) * (1 - eta);
(1 + xi) * (1 - eta);
(1 + xi) * (1 + eta);
(1 - xi) * (1 + eta)];
dN_dxi = 1/4 * [-(1 - eta), 1 - eta, 1 + eta, -(1 + eta)];
dN_deta = 1/4 * [-(1 - xi), -(1 + xi), 1 + xi, 1 - xi];
% Jacobian matrix
J = 10^-3 *[dN_dxi; dN_deta] * [0, -100; 100, -100; 100, 0; 0, 0]; % x1, y1, etc. are the nodal coordinates
% Determinant of the Jacobian matrix
detJ = det(J);
% Inverse of the Jacobian matrix
invJ = inv(J);
% B matrix (strain-displacement matrix)
B = zeros(3, 8);
for k = 1:4
B(1, 2*k-1) = dN_dxi(k);
B(2, 2*k) = dN_deta(k);
B(3, 2*k-1) = dN_deta(k);
B(3, 2*k) = dN_dxi(k);
end
% Element stiffness matrix contribution
C = E / (1 - nu^2) * [1 nu 0; nu 1 0; 0 0 (1 - nu) / 2];
Ke = t * B' * C * B * detJ;
% Add the contribution to the respective element stiffness matrix
if i == 1 && j == 1
Ke1 = Ke1 + wVec(i) * wVec(j) * Ke;
elseif i == 2 && j == 1
Ke2 = Ke2 + wVec(i) * wVec(j) * Ke;
elseif i == 2 && j == 2
Ke3 = Ke3 + wVec(i) * wVec(j) * Ke;
elseif i == 1 && j == 2
Ke4 = Ke4 + wVec(i) * wVec(j) * Ke;
end
end
end
% Display the element stiffness matrices
disp('Element Stiffness Matrix for Element 1:');
Element Stiffness Matrix for Element 1:
disp(Ke1);
1.0e+06 * 2.4027 1.1834 -1.6297 -0.4559 -0.6438 -0.3171 -0.1292 -0.4104 1.1834 2.4027 -0.4104 -0.1292 -0.3171 -0.6438 -0.4559 -1.6297 -1.6297 -0.4104 1.8368 -0.3171 0.4367 0.1100 -0.6438 0.6175 -0.4559 -0.1292 -0.3171 0.7384 0.1222 0.0346 0.6508 -0.6438 -0.6438 -0.3171 0.4367 0.1222 0.1725 0.0850 0.0346 0.1100 -0.3171 -0.6438 0.1100 0.0346 0.0850 0.1725 0.1222 0.4367 -0.1292 -0.4559 -0.6438 0.6508 0.0346 0.1222 0.7384 -0.3171 -0.4104 -1.6297 0.6175 -0.6438 0.1100 0.4367 -0.3171 1.8368
disp('Element Stiffness Matrix for Element 2:');
Element Stiffness Matrix for Element 2:
disp(Ke2);
1.0e+06 * 1.8368 0.3171 -1.6297 0.4104 -0.6438 -0.6175 0.4367 -0.1100 0.3171 0.7384 0.4559 -0.1292 -0.6508 -0.6438 -0.1222 0.0346 -1.6297 0.4559 2.4027 -1.1834 -0.1292 0.4104 -0.6438 0.3171 0.4104 -0.1292 -1.1834 2.4027 0.4559 -1.6297 0.3171 -0.6438 -0.6438 -0.6508 -0.1292 0.4559 0.7384 0.3171 0.0346 -0.1222 -0.6175 -0.6438 0.4104 -1.6297 0.3171 1.8368 -0.1100 0.4367 0.4367 -0.1222 -0.6438 0.3171 0.0346 -0.1100 0.1725 -0.0850 -0.1100 0.0346 0.3171 -0.6438 -0.1222 0.4367 -0.0850 0.1725
disp('Element Stiffness Matrix for Element 3:');
Element Stiffness Matrix for Element 3:
disp(Ke3);
1.0e+06 * 0.1725 0.0850 0.0346 0.1100 -0.6438 -0.3171 0.4367 0.1222 0.0850 0.1725 0.1222 0.4367 -0.3171 -0.6438 0.1100 0.0346 0.0346 0.1222 0.7384 -0.3171 -0.1292 -0.4559 -0.6438 0.6508 0.1100 0.4367 -0.3171 1.8368 -0.4104 -1.6297 0.6175 -0.6438 -0.6438 -0.3171 -0.1292 -0.4104 2.4027 1.1834 -1.6297 -0.4559 -0.3171 -0.6438 -0.4559 -1.6297 1.1834 2.4027 -0.4104 -0.1292 0.4367 0.1100 -0.6438 0.6175 -1.6297 -0.4104 1.8368 -0.3171 0.1222 0.0346 0.6508 -0.6438 -0.4559 -0.1292 -0.3171 0.7384
disp('Element Stiffness Matrix for Element 4:');
Element Stiffness Matrix for Element 4:
disp(Ke4);
1.0e+06 * 0.7384 0.3171 0.0346 -0.1222 -0.6438 -0.6508 -0.1292 0.4559 0.3171 1.8368 -0.1100 0.4367 -0.6175 -0.6438 0.4104 -1.6297 0.0346 -0.1100 0.1725 -0.0850 0.4367 -0.1222 -0.6438 0.3171 -0.1222 0.4367 -0.0850 0.1725 -0.1100 0.0346 0.3171 -0.6438 -0.6438 -0.6175 0.4367 -0.1100 1.8368 0.3171 -1.6297 0.4104 -0.6508 -0.6438 -0.1222 0.0346 0.3171 0.7384 0.4559 -0.1292 -0.1292 0.4104 -0.6438 0.3171 -1.6297 0.4559 2.4027 -1.1834 0.4559 -1.6297 0.3171 -0.6438 0.4104 -0.1292 -1.1834 2.4027
% (b) AssembleThe code for assembling the global stiffness matrix is as follows:
% Total number of nodes and elements
numNodes = 9; % Assuming 9 nodes in total
numElements = 4; % Number of elements
% Degree of freedom (DOF) per node
DOFperNode = 2;
% Total number of DOFs
numDOFs = numNodes * DOFperNode;
% Initialize the global stiffness matrix
KGlob = zeros(18);
% Element connectivity matrix
Elconn = [1 2 5 4; 2 3 6 5; 4 5 8 7; 5 6 9 8]; % Assuming the nodes are numbered in a specific order
% Assemble the global stiffness matrix
for element = 1:numElements
% Get the element stiffness matrix for the current element
if element == 1
Ke = Ke1;
elseif element == 2
Ke = Ke2;
elseif element == 3
Ke = Ke3;
else
Ke = Ke4;
end
% Loop over the element nodes
for i = 1:4
% Get the global node number for the current element node
node = Elconn(element, i);
% Loop over the degrees of freedom (DOFs) of the current node
for j = 1:DOFperNode
% Get the global DOF number
DOF = (node - 1) * DOFperNode + j
% Add the contribution to the global stiffness matrix
for k = 1:DOFperNode
DOF_global = (node - 1) * DOFperNode + k;
KGlob(DOF, DOF_global) = KGlob(DOF, DOF_global) + Ke(DOF, DOF_global);
end
end
end
end
DOF = 1
DOF = 2
DOF = 3
DOF = 4
DOF = 9
Index in position 1 exceeds array bounds. Index must not exceed 8.
% Display the global stiffness matrix
disp('Global Stiffness Matrix:');
disp(KGlob);
% (c) Initialize the element force vectors
fe1 = zeros(8, 1);
fe2 = zeros(8, 1);
fe3 = zeros(8, 1);
fe4 = zeros(8, 1);
% Loop over the integration points
for i = 1:dof
for j = 1:dof
xi = xiVec(i);
eta = etaVec(j);
% Shape functions and derivatives
N = 1/4 * [(1 - xi) * (1 - eta);
(1 + xi) * (1 - eta);
(1 + xi) * (1 + eta);
(1 - xi) * (1 + eta)];
% Jacobian matrix
J = [dN_dxi; dN_deta] * [x1, y1; x2, y2; x3, y3; x4, y4]; % x1, y1, etc. are the nodal coordinates
% Determinant of the Jacobian matrix
detJ = det(J);
% Element force vector contribution
fe = P * N * detJ;
% Add the contribution to the respective element force vector
if i == 1 && j == 1
fe1 = fe1 + wVec(i) * wVec(j) * fe;
elseif i == 2 && j == 1
fe2 = fe2 + wVec(i) * wVec(j) * fe;
elseif i == 2 && j == 2
fe3 = fe3 + wVec(i) * wVec(j) * fe;
elseif i == 1 && j == 2
fe4 = fe4 + wVec(i) * wVec(j) * fe;
end
end
end
% Display the element force vectors
disp('Element Force Vector for Element 1:');
disp(fe1);
disp('Element Force Vector for Element 2:');
disp(fe2);
disp('Element Force Vector for Element 3:');
disp(fe3);
disp('Element Force Vector for Element 4:');
disp(fe4);
% (D) Initialize the global force vector
F_global = zeros(numDOFs, 1);
% Assemble the global force vector
for element = 1:numElements
% Get the element force vector for the current element
if element == 1
fe = fe1;
elseif element == 2
fe = fe2;
elseif element == 3
fe = fe3;
else
fe = fe4;
end
% Loop over the element nodes
for i = 1:4
% Get the global node number for the current element node
node = Elconn(element, i);
% Loop over the degrees of freedom (DOFs) of the current node
for j = 1:DOFperNode
% Get the global DOF number
DOF = (node - 1) * DOFperNode + j;
% Add the contribution to the global force vector
F_global(DOF) = F_global(DOF) + fe(DOF);
end
end
end
% Display the global force vector
disp('Global Force Vector:');
disp(F_global);
% (e) Apply boundary conditions
% Assuming node 1 is fixed in both x and y directions
fixedDOFs = [1, 2]; % DOFs corresponding to node 1
% Set the corresponding rows and columns in the global stiffness matrix to zero, except for the diagonal entries
KGlob(fixedDOFs, :) = 0;
KGlob(:, fixedDOFs) = 0;
KGlob(fixedDOFs, fixedDOFs) = eye(length(fixedDOFs));
% Set the corresponding entries in the global force vector to zero
F_global(fixedDOFs) = 0;
% (f) Solve for the unknown displacements
U = KGlob \ F_global;
% Display the unknown displacements
disp('Unknown Displacements:');
disp(U);
% (g) Calculate the plane stress components in the center of each element
stress_center = zeros(numElements, 3); % [sigma_xx, sigma_yy, tau_xy]
for element = 1:numElements
% Get the element stiffness matrix for the current element
if element == 1
Ke = Ke1;
elseif element == 2
Ke = Ke2;
elseif element == 3
Ke = Ke3;
else
Ke = Ke4;
end
% Get the element displacements for the current element
Ue = U(Elconn(element, :));
% Calculate the strains using the strain-displacement relation
strains = B * Ue;
% Calculate the stress components using Hooke's law
stresses = C * strains;
% Store the stress components in the center of the element
stress_center(element, :) = stresses(:, 1)';
end
% Display theplane stress components in the center of each element:
disp('Plane Stress Components in the Center of Each Element:');
disp(stress_center);
  1 个评论
Dyuman Joshi
Dyuman Joshi 2023-12-26
编辑：Dyuman Joshi 2023-12-26
See the above edit, you should be able to spot where the error occurs, which is due to the value of variable DOF.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

James Tursa
James Tursa 2023-12-26
编辑：James Tursa 2023-12-26
Here is how to debug these types of errors:
Type the following at the command prompt:
dbstop if error
Then run your code. When the error occurs, the program will pause with all variables intact. Examine the variables on the line with the error to determine which ones are causing the error. Then backtrack in your code to figure out why they are not the values and/or sizes you expected.

类别

在 Help Center 和 File Exchange 中查找有关 Stress and Strain 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by