wrong matrix - provides 3x3 instead of 3x1

7 ビュー (過去 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 件のコメント
Clayton Gotberg
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 件のコメント
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 件のコメント
Cris LaPierre
Cris LaPierre 2021 年 4 月 17 日
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')

サインインしてコメントする。

Community Treasure Hunt

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

Start Hunting!

Translated by