plotting is giving me multiple lines
28 views (last 30 days)
Show older comments
Greetings,
I made some modifications to my code and I am trying to plot this, but it is giving me different lines. My intent is to plot the first row of "stress1f(:,k)" for the x axis and y axis z_1=[z(k) z(k+1)]; I tried to do it inside and outside the for loop, but matlab seems to graph the same points multiple times when i do it inside the for loop and outside it seems the same. any suggestions would be greatly appreciated it
My graph should look as the three graphs on the right when I graph the first, second, and third row against z_1=[z(k) z(k+1)];
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];
%%Z-Values
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];
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
b=stress1f(1,:)';
t=stress2f(1,:)';
tt=[b t];
z_1=[z(k) z(k+1)];
% figure(1)
% stressx=[stress1f(1,k) stress2f(1,k)];
% z_1=[z(k) z(k+1)];
% plot(stressx,z_1)
% xline(0)
% hold on
% figure(2)
% stressy=[stress1f(2,k) stress2f(2,k)];
% z_2=[z(k) z(k+1)];
% plot(stressy,z_2)
% xline(0)
% figure(3)
% stressxy=[stress1f(3,k) stress2f(3,k)];
% z_3=[z(k) z(k+1)];
% plot(stressxy,z_3)
% xline(0)
% hold on
end
figure(1)
plot(tt,z_1)
0 Comments
Answers (1)
Alan Stevens
on 19 Apr 2021
A little more like this perhaps:
%%%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];
%%Z-Values
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 Extension Coupling
shearf = A\N;
e_xf = shearf(1); e_yf = shearf(2); gamma_xyf = shearf(3);
%%Bending
bendingf = D\M;
k_xf = bendingf(1); k_yf = bendingf(2); k_xyf = bendingf(3);
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
b=stress1f(1,:)';
t=stress2f(1,:)';
tt=[b t];
z_1(k,:)=[z(k) z(k+1)];
end
z_1 = z_1';
z_1 = z_1(:);
figure(1)
stressx=[stress1f(1,:)' stress2f(1,:)']';
stressx = stressx(:);
plot(stressx,z_1)
xline(0);
figure(2)
stressy=[stress1f(2,:)' stress2f(2,:)']';
stressy = stressy(:);
plot(stressy,z_1)
xline(0);
figure(3)
stressxy=[stress1f(3,:)' stress2f(3,:)']';
stressxy = stressxy(:);
plot(stressxy,z_1)
xline(0);
0 Comments
See Also
Categories
Find more on Directed Graphs in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!