MATLAB Answers

wrong matrix - provides 3x3 instead of 3x1

3 views (last 30 days)
Eddy Ramirez
Eddy Ramirez on 17 Apr 2021
Commented: Cris LaPierre on 17 Apr 2021
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

Answers (3)

Image Analyst
Image Analyst on 17 Apr 2021
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 Comments
Clayton Gotberg
Clayton Gotberg on 17 Apr 2021
The difference between this code and the one you have posted is element-wise multiplication (A.*B) instead of matrix multiplication (A*B).

Sign in to comment.


Clayton Gotberg
Clayton Gotberg on 17 Apr 2021
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 Comments
Clayton Gotberg
Clayton Gotberg on 17 Apr 2021
Please create a new question for this and please remember to accept one of the answers if you feel it has solved your question.

Sign in to comment.


Cris LaPierre
Cris LaPierre on 17 Apr 2021
Edited: Cris LaPierre on 17 Apr 2021
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 Comments
Cris LaPierre
Cris LaPierre on 17 Apr 2021
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')

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!

Translated by