Error using symengine Unable to convert expression containing symbolic variables into double array. Apply 'subs' function first to substitute values for

6 Ansichten (letzte 30 Tage)
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
element= 1
Unable to perform assignment because value of type 'sym' is not convertible to 'double'.

Caused by:
Error using symengine
Unable to convert expression containing symbolic variables into double array. Apply 'subs' function first to substitute values for variables.
% 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

Akzeptierte Antwort

Walter Roberson
Walter Roberson am 18 Mär. 2024
Verschoben: Walter Roberson am 18 Mär. 2024
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))
kB_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
element= 1
Name Size Bytes Class Attributes A 1x1 8 double A1 1x1 8 double Ay 1x1 8 sym Ay1 1x1 8 double E 1x1 8 double Force 29x1 232 double G 1x1 8 double HAB 3x3 72 double HQB 3x3 72 double HQB_T 3x3 72 double HxB 3x3 72 double HxB_T 3x3 72 double HxB_Tc 3x3 72 double HxBc 3x3 72 double Iz 1x1 8 sym Iz1 1x1 8 double KG 6x6 288 double KKG 27x27 5832 double L 1x1 8 sym L1 1x1 8 double TAB 3x3 72 double TQB 3x3 8 sym TQB_T 3x3 8 sym alpha 1x1 8 sym alphaA 1x9 72 double alphaB 1x9 72 double angle 1x1 8 double c 1x1 8 sym c1 1x1 8 double cA 1x1 8 double cB 1x1 8 double dX 1x9 72 double dY 1x9 72 double fBA_column 3x3 8 sym fBA_curved 3x3 8 sym fBA_straight 3x3 8 sym fuq 3x3 8 sym fuq_column 3x3 8 sym fuq_straight 3x3 8 sym fux 3x3 8 sym fux_column 3x3 8 sym fux_straight 3x3 8 sym ii 1x1 8 double kB_column 3x3 72 double kB_curved 3x3 8 sym kB_straight 3x3 8 sym lA 3x3 72 double lB 3x3 72 double lambda_QB 3x3 8 sym lambda_QB_T 3x3 8 sym ndofn 1x1 8 double nnodes 1x1 8 double nu 1x1 8 double r 1x1 8 double s 1x1 8 sym s1 1x1 8 double sA 1x1 8 double sB 1x1 8 double t 2x8 128 double theta 1x1 8 double x 1x1 8 sym
Unable to perform assignment because value of type 'sym' is not convertible to 'double'.

Caused by:
Error using symengine
Unable to convert expression containing symbolic variables into double array. Apply 'subs' function first to substitute values for variables.
% 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 Kommentare
Vruddhi
Vruddhi am 23 Mär. 2024
How should i correct the above I tryed understand your comment but couldnt implement and make changes in my code. Please help
Walter Roberson
Walter Roberson am 23 Mär. 2024
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 ?

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Produkte

Community Treasure Hunt

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

Start Hunting!

Translated by