Error using symengine Unable to convert expression containing symbolic variables into double array. Apply 'subs' function first to substitute values for
조회 수: 5 (최근 30일)
이전 댓글 표시
Error in Assifinal_Vruddhi (line 177)
KG(1:3,1:3) = lA*HAB*kB_straight*HAB'*lA'; % KAA
My code is following
axes
clear
clc
close all
format short g
ndofn = 3; % number of degrees of freedom per node
% each number corresponds to a force for the corresponding DOF
Force = [0 0 0 0 0 0 0 0 0 0 -1000 0 0 -1000 0 0 -1000 0 0-1000 0 0 0 0 0 0 0 0 0 0]';
% Angles between local and global reference systems
% alphaA contains angles at the first node A
alphaA = [pi/2 pi/2 0 0 0 0 0 0 -pi/2];
% alphaB contains angles at the first node B
alphaB = [pi/2 0 0 0 0 0 0 -pi/2 -pi/2];
% distances along X and Y directions between origins of local reference systems
% in A and B; to be used in the tranlsation matrix
dX = [240 240 240 240 240 240 240 240 240]; %in
dY = [0 -240 0 0 0 0 0 -240 0]; %in
% Connectivity matrix: each column i shows the nodes that element i is connected to;
% For example, element 1 is connected to the structures nodes 1 and 2;
% In other words A==1 B==2 for the first element which is represented by
% the first column
t = [1 2 3 4 5 7 8 9;
2 3 4 5 6 8 9 10];
%for straight
syms L x A E G Ay Iz real
%translation matrix for 2D element -> simplification of 3D case
HxB = [1 0 0;
0 1 0;
0 240 1];
HxB_T = HxB'; % transpose
% applies to linear portions of bridge
% fuq is cross sectional unit flexibility matrix
fuq_straight = [1/(A*E) 0 0;
0 1/(G*Ay) 0;
0 0 1/(E*Iz)];
fux_straight = HxB_T*fuq_straight*HxB;
fBA_straight = int(fux_straight, x, 0, L);
% Numerical values given and from SAP
E = 3E7; %psi
nu = 0.29;
G = E/(2*(1+nu)); %psi
A1 = 254.469; %in^2
Iz1 = 1936; %in^4
Ay1 = 254.469; %in^2
L1 = 240; %in
% Stiffness matrix of cantilever straight element
kB_straight = eval(inv(fBA_straight));
%for columns
syms L1 x A1 E G Ay1 Iz1 real
%translation matrix for 2D element -> simplification of 3D case
HxBc = [1 0 0;
0 1 0;
0 240 1];
HxB_Tc = HxBc'; % transpose
% applies to columns
% fuq is cross sectional unit flexibility matrix
fuq_column = [1/(A1*E) 0 0;
0 1/(G*Ay1) 0;
0 0 1/(E*Iz1)];
fux_column = HxB_Tc*fuq_column*HxBc;
fBA_column = int(fux_column, x, 0, L1);
% Numerical values given and from SAP
E = 3E7; %psi
nu = 0.29;
G = E/(2*(1+nu)); %psi
A1 = 254.469; %in^2
Iz1 = 1936; %in^4
Ay1 = 254.469; %in^2
L1 = 240; %in
% Stiffness matrix of cantilever straight element
kB_column = eval(inv(fBA_column));
%for curved
syms alpha theta r E G A Ay Iz real
c = cos(alpha);
s = sin(alpha);
% rotation operator
lambda_QB_T = [c s 0;
-s c 0;
0 0 1];
lambda_QB = lambda_QB_T';
% translation operator accounts also for the different (rotated) reference
% systems at different cross sections in the curved member
HQB = [1 0 0;
0 1 0;
-240 240 1];
HQB_T = HQB';
% unitary flexibility matrix
fuq = [1/(E*A) 0 0;
0 1/(G*Ay) 0;
0 0 1/(E*Iz)];
TQB = lambda_QB_T*HQB; % Translation and rotation operator
TQB_T = TQB'; % Transpose
fux = TQB_T * fuq * TQB.*r;
fBA_curved = int(fux, alpha, 0, theta); %Flexibility matrix of curved member
theta = pi/2;
E = 3E7; %psi
nu = 0.29;
G = E/(2*(1+nu)); %psi
A = 254.469; %in^2
Iz1 = 1936; %in^4
Ay1 = 254.469; %in^2
r = 240; %in
c1 = cos(theta);
s1 = sin(theta);
% Stiffness matrix of cantilever curved element
kB_curved = eval(inv(fBA_curved));
%global matrice
nnodes = length(t)+1;
% Initializing the global stiffness matrix of the entire structure
% Its dimension is square with nnodes*ndofn rows and columns
KKG = zeros(nnodes*ndofn,nnodes*ndofn);
% Initializing the global stiffness matrix of the one element
% Its dimension is square with 2*ndofn rows and columns
KG = zeros(2*ndofn,2*ndofn);
%Initialize Matrices
HAB = zeros(3,3);
lA = zeros(3,3);
lB = zeros(3,3);
for ii=1:length(t)
HAB = [1 0 0;
0 1 0;
-dY(ii) dX(ii) 1]; % translation operator for straight members
cA = cos(alphaA(ii));
sA = sin(alphaA(ii));
cB = cos(alphaB(ii));
sB = sin(alphaB(ii));
lA = [cA -sA 0;
sA cA 0;
0 0 1];
lB = [cB -sB 0;
sB cB 0;
0 0 1];
angle = pi/2;
TAB = [ cos(angle) sin(angle) 0; % rotation/translation operator for curved members
-sin(angle) cos(angle) 0 % local ref systems in A and B are off by angle
-r*(1-cos(angle)) r*sin(angle) 1];
disp(['element= ',num2str(ii)])
disp([' '])
if ii==1 || ii==6 % for straight components
KG(1:3,1:3) = lA*HAB*kB_straight*HAB'*lA'; % KAA
KG(4:6,1:3) = -lB*kB_straight*HAB'*lA'; % KBA
KG(1:3,4:6) = -lA*HAB*kB_straight*lB'; % KAB
KG(4:6,4:6) = lB*kB_straight*lB'; % KBB
elseif ii==2 || ii==5 % for all columns
KG(1:3,1:3) = lA*HAB*kB_column*HAB'*lA'; % KAA
KG(4:6,1:3) = -lB*kB_column*HAB'*lA'; % KBA
KG(1:3,4:6) = -lA*HAB*kB_column*lB'; % KAB
KG(4:6,4:6) = lB*kB_column*lB'; % KBB
else % for curved members
KG(1:3,1:3) = lA*TAB*kB_curved*TAB'*lA'; % KAA
KG(4:6,1:3) = -lB*kB_curved*TAB'*lA'; % KBA
KG(1:3,4:6) = -lA*TAB*kB_curved*lB'; % KAB
KG(4:6,4:6) = lB*kB_curved*lB'; % KBB
end
node1 = t(1,ii);
node2 = t(2,ii);
AA = node1*3-2 : node1*3; % AA contains the dofs for the node A of the member
BB = node2*3-2 : node2*3; % BB contains the dofs for the node B of the member
% Adding the element stiffness matrix coefficients in the global
% stiffness matrix of the entire structure
KKG(AA,AA) = KKG(AA,AA) + KG(1:3,1:3);
KKG(BB,AA) = KKG(BB,AA) + KG(4:6,1:3);
KKG(AA,BB) = KKG(AA,BB) + KG(1:3,4:6);
KKG(BB,BB) = KKG(BB,BB) + KG(4:6,4:6);
end
% Eliminate rows and columns corresponding to constrained DOFs
Force([1:3,28:30],:) = [];
KK = KKG([1,3:6,10:12,16:19,21],[1,3:6,10:12,16:19,21]);
displ = KK\Force;
%Displacements of unconstrained nodes
ux = displ([4,7,10,13,16,19,22,25]);
uy = displ([11,14,17,20]);
phi = displ([6,9,12,15,18,21,24]);
% Table of nodal displacements for the structure
disp([' Node ' ' ux [in] ' 'uy [in] ' 'phi [rad]'])
disp([[2;4;6;7] , ux , [0;uy;0] , phi])
% Displacement components for all nodes
displALL_x = [ux(1),ux(2),0,ux(3),0,ux(4),ux(5)];
displALL_y = [0,uy(1),0,uy(2),0,uy(3),0];
n_coord = [000 120 120 195 270 270 390;...
035 075 000 085 000 075 035].*12;
xcoord = n_coord(1,:);
ycoord = n_coord(2,:);
scale = 100;
figure(1), plot(xcoord,ycoord,'ks')
hold on
plot(xcoord + scale.*displALL_x , ycoord + scale.*displALL_y,'ro')
legend('undeformed nodes','deformed nodes')
xlabel('X axis [in]')
ylabel('Y axis [in]')
axis equal
grid on
댓글 수: 0
채택된 답변
Walter Roberson
2024년 3월 18일
이동: Walter Roberson
2024년 3월 18일
axes
clear
clc
close all
format short g
ndofn = 3; % number of degrees of freedom per node
% each number corresponds to a force for the corresponding DOF
Force = [0 0 0 0 0 0 0 0 0 0 -1000 0 0 -1000 0 0 -1000 0 0-1000 0 0 0 0 0 0 0 0 0 0]';
% Angles between local and global reference systems
% alphaA contains angles at the first node A
alphaA = [pi/2 pi/2 0 0 0 0 0 0 -pi/2];
% alphaB contains angles at the first node B
alphaB = [pi/2 0 0 0 0 0 0 -pi/2 -pi/2];
% distances along X and Y directions between origins of local reference systems
% in A and B; to be used in the tranlsation matrix
dX = [240 240 240 240 240 240 240 240 240]; %in
dY = [0 -240 0 0 0 0 0 -240 0]; %in
% Connectivity matrix: each column i shows the nodes that element i is connected to;
% For example, element 1 is connected to the structures nodes 1 and 2;
% In other words A==1 B==2 for the first element which is represented by
% the first column
t = [1 2 3 4 5 7 8 9;
2 3 4 5 6 8 9 10];
%for straight
syms L x A E G Ay Iz real
%translation matrix for 2D element -> simplification of 3D case
HxB = [1 0 0;
0 1 0;
0 240 1];
HxB_T = HxB'; % transpose
% applies to linear portions of bridge
% fuq is cross sectional unit flexibility matrix
fuq_straight = [1/(A*E) 0 0;
0 1/(G*Ay) 0;
0 0 1/(E*Iz)];
fux_straight = HxB_T*fuq_straight*HxB;
fBA_straight = int(fux_straight, x, 0, L);
fBA_straight is the integral to the unbound symbolic variable L, so fBA_straight is going to be symbolic.
% Numerical values given and from SAP
E = 3E7; %psi
nu = 0.29;
G = E/(2*(1+nu)); %psi
A1 = 254.469; %in^2
Iz1 = 1936; %in^4
Ay1 = 254.469; %in^2
L1 = 240; %in
% Stiffness matrix of cantilever straight element
kB_straight = eval(inv(fBA_straight))
You assigned to A1 and L1 and ay1 not to A and L and ay, so kB_straight remains symbolic after the eval() step.
%for columns
syms L1 x A1 E G Ay1 Iz1 real
%translation matrix for 2D element -> simplification of 3D case
HxBc = [1 0 0;
0 1 0;
0 240 1];
HxB_Tc = HxBc'; % transpose
% applies to columns
% fuq is cross sectional unit flexibility matrix
fuq_column = [1/(A1*E) 0 0;
0 1/(G*Ay1) 0;
0 0 1/(E*Iz1)];
fux_column = HxB_Tc*fuq_column*HxBc;
fBA_column = int(fux_column, x, 0, L1);
% Numerical values given and from SAP
E = 3E7; %psi
nu = 0.29;
G = E/(2*(1+nu)); %psi
A1 = 254.469; %in^2
Iz1 = 1936; %in^4
Ay1 = 254.469; %in^2
L1 = 240; %in
% Stiffness matrix of cantilever straight element
kB_column = eval(inv(fBA_column));
%for curved
syms alpha theta r E G A Ay Iz real
c = cos(alpha);
s = sin(alpha);
% rotation operator
lambda_QB_T = [c s 0;
-s c 0;
0 0 1];
lambda_QB = lambda_QB_T';
% translation operator accounts also for the different (rotated) reference
% systems at different cross sections in the curved member
HQB = [1 0 0;
0 1 0;
-240 240 1];
HQB_T = HQB';
% unitary flexibility matrix
fuq = [1/(E*A) 0 0;
0 1/(G*Ay) 0;
0 0 1/(E*Iz)];
TQB = lambda_QB_T*HQB; % Translation and rotation operator
TQB_T = TQB'; % Transpose
fux = TQB_T * fuq * TQB.*r;
fBA_curved = int(fux, alpha, 0, theta); %Flexibility matrix of curved member
theta = pi/2;
E = 3E7; %psi
nu = 0.29;
G = E/(2*(1+nu)); %psi
A = 254.469; %in^2
Iz1 = 1936; %in^4
Ay1 = 254.469; %in^2
r = 240; %in
c1 = cos(theta);
s1 = sin(theta);
% Stiffness matrix of cantilever curved element
kB_curved = eval(inv(fBA_curved));
%global matrice
nnodes = length(t)+1;
% Initializing the global stiffness matrix of the entire structure
% Its dimension is square with nnodes*ndofn rows and columns
KKG = zeros(nnodes*ndofn,nnodes*ndofn);
% Initializing the global stiffness matrix of the one element
% Its dimension is square with 2*ndofn rows and columns
KG = zeros(2*ndofn,2*ndofn);
You assign KG as numeric zeros.
%Initialize Matrices
HAB = zeros(3,3);
lA = zeros(3,3);
lB = zeros(3,3);
for ii=1:length(t)
HAB = [1 0 0;
0 1 0;
-dY(ii) dX(ii) 1]; % translation operator for straight members
cA = cos(alphaA(ii));
sA = sin(alphaA(ii));
cB = cos(alphaB(ii));
sB = sin(alphaB(ii));
lA = [cA -sA 0;
sA cA 0;
0 0 1];
lB = [cB -sB 0;
sB cB 0;
0 0 1];
angle = pi/2;
TAB = [ cos(angle) sin(angle) 0; % rotation/translation operator for curved members
-sin(angle) cos(angle) 0 % local ref systems in A and B are off by angle
-r*(1-cos(angle)) r*sin(angle) 1];
disp(['element= ',num2str(ii)])
disp([' '])
whos
if ii==1 || ii==6 % for straight components
KG(1:3,1:3) = lA*HAB*kB_straight*HAB'*lA'; % KAA
You use the symbolic kB_straight on the right-hand side, so the right hand side comes out symbolic involving unbound variables. But the left-hand side is purely numeric, and numeric arrays cannot hold expressions involving unbound symbolic variables.
KG(4:6,1:3) = -lB*kB_straight*HAB'*lA'; % KBA
KG(1:3,4:6) = -lA*HAB*kB_straight*lB'; % KAB
KG(4:6,4:6) = lB*kB_straight*lB'; % KBB
elseif ii==2 || ii==5 % for all columns
KG(1:3,1:3) = lA*HAB*kB_column*HAB'*lA'; % KAA
KG(4:6,1:3) = -lB*kB_column*HAB'*lA'; % KBA
KG(1:3,4:6) = -lA*HAB*kB_column*lB'; % KAB
KG(4:6,4:6) = lB*kB_column*lB'; % KBB
else % for curved members
KG(1:3,1:3) = lA*TAB*kB_curved*TAB'*lA'; % KAA
KG(4:6,1:3) = -lB*kB_curved*TAB'*lA'; % KBA
KG(1:3,4:6) = -lA*TAB*kB_curved*lB'; % KAB
KG(4:6,4:6) = lB*kB_curved*lB'; % KBB
end
node1 = t(1,ii);
node2 = t(2,ii);
AA = node1*3-2 : node1*3; % AA contains the dofs for the node A of the member
BB = node2*3-2 : node2*3; % BB contains the dofs for the node B of the member
% Adding the element stiffness matrix coefficients in the global
% stiffness matrix of the entire structure
KKG(AA,AA) = KKG(AA,AA) + KG(1:3,1:3);
KKG(BB,AA) = KKG(BB,AA) + KG(4:6,1:3);
KKG(AA,BB) = KKG(AA,BB) + KG(1:3,4:6);
KKG(BB,BB) = KKG(BB,BB) + KG(4:6,4:6);
end
% Eliminate rows and columns corresponding to constrained DOFs
Force([1:3,28:30],:) = [];
KK = KKG([1,3:6,10:12,16:19,21],[1,3:6,10:12,16:19,21]);
displ = KK\Force;
%Displacements of unconstrained nodes
ux = displ([4,7,10,13,16,19,22,25]);
uy = displ([11,14,17,20]);
phi = displ([6,9,12,15,18,21,24]);
% Table of nodal displacements for the structure
disp([' Node ' ' ux [in] ' 'uy [in] ' 'phi [rad]'])
disp([[2;4;6;7] , ux , [0;uy;0] , phi])
% Displacement components for all nodes
displALL_x = [ux(1),ux(2),0,ux(3),0,ux(4),ux(5)];
displALL_y = [0,uy(1),0,uy(2),0,uy(3),0];
n_coord = [000 120 120 195 270 270 390;...
035 075 000 085 000 075 035].*12;
xcoord = n_coord(1,:);
ycoord = n_coord(2,:);
scale = 100;
figure(1), plot(xcoord,ycoord,'ks')
hold on
plot(xcoord + scale.*displALL_x , ycoord + scale.*displALL_y,'ro')
legend('undeformed nodes','deformed nodes')
xlabel('X axis [in]')
ylabel('Y axis [in]')
axis equal
grid on
댓글 수: 8
Walter Roberson
2024년 3월 23일
You have
ux = displ([2,3,8,9]);
which makes ux a vector of length 4.
You tried to do
disp([[2;3;4;5;6;7;8] , ux , [0;uy;0] , phi])
which would require that ux is a vector of length 7.
You tried to do
displALL_x = [ux(2),ux(3),uy(7),ux(8)];
which would have required that ux is a vector of length at least 8.
These definitions of ux are incompatible.
You also have
uy = displ([2,3,4,5,6,7,8]);
which makes uy have a vector of length 7. It is questionable whether it makes sense to extract the same elements of displ for both ux and uy.
Getting back to
disp([[2;3;4;5;6;7;8] , ux , [0;uy;0] , phi])
this requires that uy is a vector of length 5, which is incompatible with the way it is initialized.
I cannot tell you how to repair these problems; you need to re-think what ux and uy mean.
Question: in
displALL_x = [ux(2),ux(3),uy(7),ux(8)];
are you sure that you want uy(7) in a vector that is intended to be ALL_x ?
추가 답변 (0개)
참고 항목
카테고리
Help Center 및 File Exchange에서 Mathematics에 대해 자세히 알아보기
제품
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!