wrong matrix - provides 3x3 instead of 3x1

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

Answer 1

[Image Analyst]

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;

[Eddy Ramirez]

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]

The difference between this code and the one you have posted is element-wise multiplication (A.*B) instead of matrix multiplication (A*B).

Answer 2

[Clayton Gotberg]

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.

[Eddy Ramirez]

that did the trick removing the . out of the formulas

thank you

[Eddy Ramirez]

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]

Please create a new question for this and please remember to accept one of the answers if you feel it has solved your question.

Answer 3

[Cris LaPierre]

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

[Eddy Ramirez]

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]

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')