wrong matrix - provides 3x3 instead of 3x1

14 次查看（过去 30 天）
Eddy Ramirez
Eddy Ramirez 2021-4-17
评论： Cris LaPierre 2021-4-17
Greetings,
I am running the code below and for the stress I am getting a 3x3 and I am not sure what the issue could be. I tried to run "sym" and created a 3x1 matrix
stress1=stress==Q_bar{k}.*shearf+z(k).*Q_bar{k}.*bendingf;%%BOTTOM LAYER
But this equation does not work either, and I cant find the reason behind it
close all
clear all
clc
%%%FIBER
Ef=220e9;%[N/m] GPA to Newton/square meter
Vf=.63; %fiber volume fraction
vf=.33;%fiber poissions ratio
%%MATRIX
Em=10e9;%[N/m] GPA to Newton/square meter
Vm=1-Vf; %%matrix volume fraction
vm=.33;%%matrix poissons ratio
Em2=Em/(1-vm^2); %equation 3.31
%%Lamina's Thickness
h=1e-3;%[N/m] mm to Newton/square meter
%%Shear modulus
Gf=Ef/(2*(1+vf));
Gm=Em/(2*(1+vm));
%%Lamina Properties
E1=(Vf*Ef)+(Vm*Em);
E2=(Ef*Em)/(Vf*Gm+Vm*Gf);
E_2=Ef*Em2/(Vf*Em2+Vm*Ef); %equation 3.32
v12=(Vf*vf)+(Vm*vm);
v21=(E2/E1)*v12;
G12=(Gf*Gm)/((Vf*Gm)+(Vm*Gf));
% Reduced local in plane stiffness Q
Q11=E1/(1-v12*v21);
Q12=(v21*E1)/(1-v12*v21);
Q22=E2/(1-v12*v21);
Q66=G12;
Q=[Q11,Q12,0;
Q12,Q22,0;
0,0,Q66];
%%Laminate rotations in degrees
stheta=[0,45,-45,90,30,60,-30,-60,-60,-30,60,30,90,-45,45,0];
z=(-length(stheta)/2+(0:length(stheta)))'*h;
% Global Q matrix:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:length(stheta)
m=cosd(stheta(i));
n=sind(stheta(i));
T=[m.^2,n.^2,2.*m.*n;
n.^2,m.^2,-2.*m.*n;
-m.*n,m.*n,m.^2-n.^2];
% Global transformed reduced stiffness coefficients Q_bar
Q__bar=inv(T).*Q.*T;
Q_bar{i}=Q__bar;
end
Aij=0;
Bij=0;
Dij=0;
for j=1:length(stheta)
A_j=Q_bar{j}*(z(j+1)-z(j));
Aj{j}=A_j;
Aij=Aij+Aj{j};
B_j=Q_bar{j}*((z(j+1))^2-(z(j))^2);
Bj{j}=B_j;
Bij=Bij+Bj{j};
D_j=Q_bar{j}*((z(j+1))^3-(z(j))^3);
Dj{j}=D_j;
Dij=Dij+Dj{j};
end
A=Aij;
B=(1/2)*Bij;
D=(1/3)*Dij;
%%Forces
Nx=2e6;
Ny=4.6e6;
Ns=0;
N=[Nx;Ny;Ns];
%%Moments
Mx=3e3;
My=0;
Ms=-1e-3;
M=[Mx;My;Ms];
%%Shear
e_x=sym('Epsilonx');
e_y=sym('Epsilony');
gamma_xy=sym('Gammaxy');
shear=[e_x;e_y; gamma_xy];
%%Bending Twist
k_x=sym('Kx');
k_y=sym('Ky');
k_xy=sym('Kxy');
bending=[k_x;k_y;k_xy];
%%Shear Extension Coupling
SEC=N==A*shear;
SEC_F=solve(SEC);
e_xf=vpa(SEC_F.Epsilonx);
e_yf=vpa(SEC_F.Epsilony);
gamma_xyf=(SEC_F.Gammaxy);
shearf=[e_xf;e_yf; gamma_xyf];
%%Bending
BEC=M==D*bending;
BEC_F=solve(BEC);
k_xf=vpa(BEC_F.Kx);
k_yf=vpa(BEC_F.Ky);
k_xyf=vpa(BEC_F.Kxy);
bendingf=[k_xf;k_yf;k_xyf];
test=Q_bar{1}.*shearf+z(1,1)*Q_bar{1}.*bendingf;
test1=Q_bar{1}.*shearf+z(2,1)*Q_bar{1}.*bendingf;
test2=Q_bar{2};
for k=1:length(stheta)
stress1=Q_bar{k}.*shearf+z(k).*Q_bar{k}.*bendingf;%%BOTTOM LAYER
stress2=Q_bar{k}.*shearf+z(k+1).*Q_bar{k}.*bendingf;%%TOP LAYER
stress1f{k}=stress1; %Bottom
stress2f{k}=stress2; %Top
end
keyboard

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

回答（3 个）

Image Analyst
Image Analyst 2021-4-17
Well isn't Q_bar a 3x3 matrix? So of course stress1 would also be 3x3.
And this is bad in terms of readability:
Q__bar=inv(T).*Q.*T; % Double underlines - hard to see that!
Q_bar{i}=Q__bar;
Simply do
Q_bar{i} = inv(T).*Q.*T;
  2 个评论
Eddy Ramirez
Eddy Ramirez 2021-4-17
hmm I have ran other codes and I am able to get a 3x1 even with having a 3x3 matrix though
clear all;
clc;
%%DATA%%
E1=150;
E2=14;
E3=E2;
G12=3.7;
G23=1.9;
G13=G23;
V12=0.33;
V13=0.37;
V23=V13;
V31=V13;
%%%3-DIMENSIONAL COORDINATES%%%
S11=1/E1;
S21=(-V12)/E1;
S12=S21;
S31=(-V13)/E1;
S13=S31;
S22=1/E2;
S32=(-V23)/E2;
S23=S32;
S33=1/E3;
S44=1/G23;
S55=1/G13;
S66=1/G12;
format shortE
S=[S11 S12 S13 0 0 0;
S12 S22 S23 0 0 0;
S13 S23 S33 0 0 0;
0 0 0 S44 0 0;
0 0 0 0 S55 0;
0 0 0 0 0 S66];
C=inv(S);
%%%TO FIND STRESS & STRAIN IN ALL PRINCIPLE DIRECTIONS (3x3 MATRIX)
sigma1=0;
sigma2=sym('sigma_2');
sigma3=sym('sigma_3');
stress=[sigma1; sigma2; sigma3];
eps1=sym('epsilon_1');
eps2=.7/100;
eps3=0;
strain=[eps1; eps2; eps3];
S_3x3=[S11,S12,S13;S21, S22, S23; S31, S23, S33];
equation=strain==S_3x3*stress;
solution=solve(equation);
Epsilonf1=vpa(solution.epsilon_1);
Sigmaf2=vpa(solution.sigma_2*1e-6);
Sigmaf3=vpa(solution.sigma_3*1e-6);
keyboard
Clayton Gotberg
Clayton Gotberg 2021-4-17
编辑：Clayton Gotberg 2021-4-17
The difference between this code and the one you have posted is element-wise multiplication (A.*B) instead of matrix multiplication (A*B).

请先登录，再进行评论。

Clayton Gotberg
Clayton Gotberg 2021-4-17
When you're multiplying matrices, you several times use element-wise multiplication instead of matrix multiplication.
A = [1 2 3; 4 5 6; 7 8 9];
B = [1;2;3];
C = A*B; % -> C is a 3x1 matrix [14;32;50];
D = A.*B; % -> D is a 3x3 matrix [1 2 3; 8 10 12; 21 24 27];
If you expect Q_bar to contain 3x3 matrices, you need to switch from element-wise multiplication.
  3 个评论
显示 1更早的评论
Eddy Ramirez
Eddy Ramirez 2021-4-17
what would be the best to plot it? I did the following, but I get a blank graph
close all
clear all
clc
%%%FIBER
Ef=220e9;%[N/m] GPA to Newton/square meter
Vf=.63; %fiber volume fraction
vf=.33;%fiber poissions ratio
%%MATRIX
Em=10e9;%[N/m] GPA to Newton/square meter
Vm=1-Vf; %%matrix volume fraction
vm=.33;%%matrix poissons ratio
Em2=Em/(1-vm^2); %equation 3.31
%%Lamina's Thickness
h=1e-3;%[N/m] mm to Newton/square meter
%%Shear modulus
Gf=Ef/(2*(1+vf));
Gm=Em/(2*(1+vm));
%%Lamina Properties
E1=(Vf*Ef)+(Vm*Em);
E2=(Ef*Em)/(Vf*Gm+Vm*Gf);
E_2=Ef*Em2/(Vf*Em2+Vm*Ef); %equation 3.32
v12=(Vf*vf)+(Vm*vm);
v21=(E2/E1)*v12;
G12=(Gf*Gm)/((Vf*Gm)+(Vm*Gf));
% Reduced local in plane stiffness Q
Q11=E1/(1-v12*v21);
Q12=(v21*E1)/(1-v12*v21);
Q22=E2/(1-v12*v21);
Q66=G12;
Q=[Q11,Q12,0;
Q12,Q22,0;
0,0,Q66];
%%Laminate rotations in degrees
theta=[0,45,-45,90,30,60,-30,-60,-60,-30,60,30,90,-45,45,0];
%%Laminate rotations in radians
% theta=pi/180*[0;45;-45;90;30;60;-30;-60]';%%DEGREES TO RADIANS
% theta=[theta;theta(end:-1:1)];
z=(-length(theta)/2+(0:length(theta)))'*h;
% Global Q matrix:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:length(theta)
m=cosd(theta(i));
n=sind(theta(i));
T=[m.^2,n.^2,2.*m.*n;
n.^2,m.^2,-2.*m.*n;
-m.*n,m.*n,m.^2-n.^2];
% Global transformed reduced stiffness coefficients Q_bar
Q_bar{i}=inv(T).*Q.*T;
end
Aij=0;
Bij=0;
Dij=0;
for j=1:length(theta)
A_j=Q_bar{j}*(z(j+1)-z(j));
Aj{j}=A_j;
Aij=Aij+Aj{j};
B_j=Q_bar{j}*((z(j+1))^2-(z(j))^2);
Bj{j}=B_j;
Bij=Bij+Bj{j};
D_j=Q_bar{j}*((z(j+1))^3-(z(j))^3);
Dj{j}=D_j;
Dij=Dij+Dj{j};
end
A=Aij;
B=(1/2)*Bij;
D=(1/3)*Dij;
%%Forces
Nx=2e6;
Ny=4.6e6;
Ns=0;
N=[Nx;Ny;Ns];
%%Moments
Mx=3e3;
My=0;
Ms=-1e-3;
M=[Mx;My;Ms];
%%Shear
e_x=sym('Epsilonx');
e_y=sym('Epsilony');
gamma_xy=sym('Gammaxy');
shear=[e_x;e_y; gamma_xy];
%%Bending Twist
k_x=sym('Kx');
k_y=sym('Ky');
k_xy=sym('Kxy');
bending=[k_x;k_y;k_xy];
%%Shear Extension Coupling
SEC=N==A*shear;
SEC_F=solve(SEC);
e_xf=vpa(SEC_F.Epsilonx);
e_yf=vpa(SEC_F.Epsilony);
gamma_xyf=(SEC_F.Gammaxy);
shearf=[e_xf;e_yf; gamma_xyf];
%%Bending
BEC=M==D*bending;
BEC_F=solve(BEC);
k_xf=vpa(BEC_F.Kx);
k_yf=vpa(BEC_F.Ky);
k_xyf=vpa(BEC_F.Kxy);
bendingf=[k_xf;k_yf;k_xyf];
test=Q_bar{1}*shearf+z(1,1)*Q_bar{1}*bendingf;
test1=Q_bar{1}.*shearf+z(2,1)*Q_bar{1}.*bendingf;
% plot(test(1,1), z(1,1));
for k=1:length(theta)
stress1=Q_bar{k}*shearf+z(k)*Q_bar{k}*bendingf;%%BOTTOM LAYER
stress2=Q_bar{k}*shearf+z(k+1)*Q_bar{k}*bendingf;%%TOP LAYER
stress1f{k}=stress1; %Bottom
stress2f{k}=stress2; %Top
end
% for p=1:length(theta)
% test=stress1f{p,1};
% end
b0=stress1f{1};
b1=stress1f{2};
b2=stress1f{3};
b3=stress1f{4};
b4=stress1f{5};
b5=stress1f{6};
b6=stress1f{1};
b7=stress1f{7};
b8=stress1f{8};
b9=stress1f{9};
b10=stress1f{10};
b11=stress1f{11};
b12=stress1f{12};
b13=stress1f{13};
b14=stress1f{14};
b15=stress1f{15};
t0=stress2f{1};
t1=stress2f{2};
t2=stress2f{3};
t3=stress2f{4};
t4=stress2f{5};
t5=stress2f{6};
t6=stress2f{1};
t7=stress2f{7};
t8=stress2f{8};
t9=stress2f{9};
t10=stress2f{10};
t11=stress2f{11};
t12=stress2f{12};
t13=stress2f{13};
t14=stress2f{14};
t15=stress2f{15};
plot(b0(1,1), z(1,1))
keyboard
hold on
plot(t0(1,1), z(2,1))
plot(b1(1,1), z(2,1))
plot(t1(1,1), z(3,1))
plot(b2(1,1), z(3,1))
plot(t2(1,1), z(4,1))
plot(b3(1,1), z(4,1))
plot(t3(1,1), z(5,1))
plot(b4(1,1), z(5,1))
plot(t4(1,1), z(6,1))
plot(b5(1,1), z(6,1))
plot(t5(1,1), z(7,1))
plot(b6(1,1), z(7,1))
plot(t6(1,1), z(8,1))
plot(b7(1,1), z(8,1))
plot(t7(1,1), z(9,1))
plot(b8(1,1), z(9,1))
plot(t8(1,1), z(10,1))
plot(b9(1,1), z(10,1))
plot(t9(1,1), z(11,1))
plot(b10(1,1), z(11,1))
plot(t10(1,1), z(12,1))
plot(b11(1,1), z(12,1))
plot(t11(1,1), z(13,1))
plot(b12(1,1), z(13,1))
plot(t12(1,1), z(14,1))
plot(b13(1,1), z(14,1))
plot(t13(1,1), z(15,1))
plot(b14(1,1), z(15,1))
plot(t14(1,1), z(16,1))
plot(b15(1,1), z(16,1))
plot(t15(1,1), z(17,1))
hold off
keyboard
Clayton Gotberg
Clayton Gotberg 2021-4-17
Please create a new question for this and please remember to accept one of the answers if you feel it has solved your question.

请先登录，再进行评论。

Cris LaPierre
Cris LaPierre 2021-4-17
编辑：Cris LaPierre 2021-4-17
Follow the dimensions of your variables to track it down.
  • Q_bar{k} is 3x3
  • Q_bar{k} = inv(T).*Q.*T where Q and T are 3x3
  • shearf is 3x1
  • bendingf is 3x1
Perhaps you don't want to be doing elementwise multiplication in your calculation of stresses. When you perform matrix multiplication, a 3x3 * 3x1 = 3x1. When you do elementwise, a 3x3 .* 3x1 = 3x3.
for k=1:length(stheta)
stress1=Q_bar{k}*shearf+z(k)*Q_bar{k}*bendingf;%%BOTTOM LAYER
% ^ ^ ^ perform matrix multiplication
stress2=Q_bar{k}*shearf+z(k+1)*Q_bar{k}*bendingf;%%TOP LAYER
% ^ ^ ^ perform matrix multiplication
stress1f{k}=stress1; %Bottom
stress2f{k}=stress2; %Top
end
  2 个评论
Eddy Ramirez
Eddy Ramirez 2021-4-17
thank you any idea how I can run the plot here I get a blank plot
D = A.*B; % -> D is a 3x3 matrix [1 2 3; 8 10 12; 21 24 27];
If you expect Q_bar to contain 3x3 matrices, you need to switch from element-wise multiplication.
2 Comments
Show 1 older comment
Eddy Ramirez
Eddy Ramirez less than a minute ago
what would be the best to plot it? I did the following, but I get a blank graph
close all
clear all
clc
%%%FIBER
Ef=220e9;%[N/m] GPA to Newton/square meter
Vf=.63; %fiber volume fraction
vf=.33;%fiber poissions ratio
%%MATRIX
Em=10e9;%[N/m] GPA to Newton/square meter
Vm=1-Vf; %%matrix volume fraction
vm=.33;%%matrix poissons ratio
Em2=Em/(1-vm^2); %equation 3.31
%%Lamina's Thickness
h=1e-3;%[N/m] mm to Newton/square meter
%%Shear modulus
Gf=Ef/(2*(1+vf));
Gm=Em/(2*(1+vm));
%%Lamina Properties
E1=(Vf*Ef)+(Vm*Em);
E2=(Ef*Em)/(Vf*Gm+Vm*Gf);
E_2=Ef*Em2/(Vf*Em2+Vm*Ef); %equation 3.32
v12=(Vf*vf)+(Vm*vm);
v21=(E2/E1)*v12;
G12=(Gf*Gm)/((Vf*Gm)+(Vm*Gf));
% Reduced local in plane stiffness Q
Q11=E1/(1-v12*v21);
Q12=(v21*E1)/(1-v12*v21);
Q22=E2/(1-v12*v21);
Q66=G12;
Q=[Q11,Q12,0;
Q12,Q22,0;
0,0,Q66];
%%Laminate rotations in degrees
theta=[0,45,-45,90,30,60,-30,-60,-60,-30,60,30,90,-45,45,0];
%%Laminate rotations in radians
% theta=pi/180*[0;45;-45;90;30;60;-30;-60]';%%DEGREES TO RADIANS
% theta=[theta;theta(end:-1:1)];
z=(-length(theta)/2+(0:length(theta)))'*h;
% Global Q matrix:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i=1:length(theta)
m=cosd(theta(i));
n=sind(theta(i));
T=[m.^2,n.^2,2.*m.*n;
n.^2,m.^2,-2.*m.*n;
-m.*n,m.*n,m.^2-n.^2];
% Global transformed reduced stiffness coefficients Q_bar
Q_bar{i}=inv(T).*Q.*T;
end
Aij=0;
Bij=0;
Dij=0;
for j=1:length(theta)
A_j=Q_bar{j}*(z(j+1)-z(j));
Aj{j}=A_j;
Aij=Aij+Aj{j};
B_j=Q_bar{j}*((z(j+1))^2-(z(j))^2);
Bj{j}=B_j;
Bij=Bij+Bj{j};
D_j=Q_bar{j}*((z(j+1))^3-(z(j))^3);
Dj{j}=D_j;
Dij=Dij+Dj{j};
end
A=Aij;
B=(1/2)*Bij;
D=(1/3)*Dij;
%%Forces
Nx=2e6;
Ny=4.6e6;
Ns=0;
N=[Nx;Ny;Ns];
%%Moments
Mx=3e3;
My=0;
Ms=-1e-3;
M=[Mx;My;Ms];
%%Shear
e_x=sym('Epsilonx');
e_y=sym('Epsilony');
gamma_xy=sym('Gammaxy');
shear=[e_x;e_y; gamma_xy];
%%Bending Twist
k_x=sym('Kx');
k_y=sym('Ky');
k_xy=sym('Kxy');
bending=[k_x;k_y;k_xy];
%%Shear Extension Coupling
SEC=N==A*shear;
SEC_F=solve(SEC);
e_xf=vpa(SEC_F.Epsilonx);
e_yf=vpa(SEC_F.Epsilony);
gamma_xyf=(SEC_F.Gammaxy);
shearf=[e_xf;e_yf; gamma_xyf];
%%Bending
BEC=M==D*bending;
BEC_F=solve(BEC);
k_xf=vpa(BEC_F.Kx);
k_yf=vpa(BEC_F.Ky);
k_xyf=vpa(BEC_F.Kxy);
bendingf=[k_xf;k_yf;k_xyf];
test=Q_bar{1}*shearf+z(1,1)*Q_bar{1}*bendingf;
test1=Q_bar{1}.*shearf+z(2,1)*Q_bar{1}.*bendingf;
% plot(test(1,1), z(1,1));
for k=1:length(theta)
stress1=Q_bar{k}*shearf+z(k)*Q_bar{k}*bendingf;%%BOTTOM LAYER
stress2=Q_bar{k}*shearf+z(k+1)*Q_bar{k}*bendingf;%%TOP LAYER
stress1f{k}=stress1; %Bottom
stress2f{k}=stress2; %Top
end
b0=stress1f{1};
b1=stress1f{2};
b2=stress1f{3};
b3=stress1f{4};
b4=stress1f{5};
b5=stress1f{6};
b6=stress1f{1};
b7=stress1f{7};
b8=stress1f{8};
b9=stress1f{9};
b10=stress1f{10};
b11=stress1f{11};
b12=stress1f{12};
b13=stress1f{13};
b14=stress1f{14};
b15=stress1f{15};
t0=stress2f{1};
t1=stress2f{2};
t2=stress2f{3};
t3=stress2f{4};
t4=stress2f{5};
t5=stress2f{6};
t6=stress2f{1};
t7=stress2f{7};
t8=stress2f{8};
t9=stress2f{9};
t10=stress2f{10};
t11=stress2f{11};
t12=stress2f{12};
t13=stress2f{13};
t14=stress2f{14};
t15=stress2f{15};
plot(b0(1,1), z(1,1))
hold on
plot(t0(1,1), z(2,1))
plot(b1(1,1), z(2,1))
plot(t1(1,1), z(3,1))
plot(b2(1,1), z(3,1))
plot(t2(1,1), z(4,1))
plot(b3(1,1), z(4,1))
plot(t3(1,1), z(5,1))
plot(b4(1,1), z(5,1))
plot(t4(1,1), z(6,1))
plot(b5(1,1), z(6,1))
plot(t5(1,1), z(7,1))
plot(b6(1,1), z(7,1))
plot(t6(1,1), z(8,1))
plot(b7(1,1), z(8,1))
plot(t7(1,1), z(9,1))
plot(b8(1,1), z(9,1))
plot(t8(1,1), z(10,1))
plot(b9(1,1), z(10,1))
plot(t9(1,1), z(11,1))
plot(b10(1,1), z(11,1))
plot(t10(1,1), z(12,1))
plot(b11(1,1), z(12,1))
plot(t11(1,1), z(13,1))
plot(b12(1,1), z(13,1))
plot(t12(1,1), z(14,1))
plot(b13(1,1), z(14,1))
plot(t13(1,1), z(15,1))
plot(b14(1,1), z(15,1))
plot(t14(1,1), z(16,1))
plot(b15(1,1), z(16,1))
plot(t15(1,1), z(17,1))
hold off
Cris LaPierre
Cris LaPierre 2021-4-17
Ran in:
When you plot a single point, you must specify a marker in order to see it. By default, MATLAB does not include one. You can see the available options here.
% no marker specified, so figure appears blank
plot(1,2)
figure
% Marker specified, so can see individual points
plot(1,2,'o')

请先登录，再进行评论。

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by