wrong matrix - provides 3x3 instead of 3x1
10 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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
0 commentaires
Réponses (3)
Image Analyst
le 17 Avr 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 commentaires
Clayton Gotberg
le 17 Avr 2021
Modifié(e) : Clayton Gotberg
le 17 Avr 2021
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
le 17 Avr 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 commentaires
Clayton Gotberg
le 17 Avr 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.
Cris LaPierre
le 17 Avr 2021
Modifié(e) : Cris LaPierre
le 17 Avr 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
Voir également
Catégories
En savoir plus sur Line Plots dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!