Issues with plotting/best method
1 visualización (últimos 30 días)
Mostrar comentarios más antiguos
Greetings,
I would like to plot the following, but I get a blank graph and I am not sure if it is the way I am inputting the values
any recommendations? or a different approach?
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))
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
8 comentarios
Clayton Gotberg
el 17 de Abr. de 2021
To connect the dots, just plot all of the points at once in the order you want them to be connected, using a line.
If I have three points (x_a, y_a), (x_b,y_b), and (x_c,y_c), I can plot them separately or together:
%Separately
plot(x_a, y_a,'k.')
plot(x_b, y_b,'k.')
plot(x_c, y_c,'k.')
%Together
plot([x_a x_b x_c],[y_a y_b y_c],'k-')
% Or
X = [x_a x_b x_c];
Y = [y_a y_b y_c];
plot(X,Y,'k-')
Respuestas (2)
LO
el 17 de Abr. de 2021
short answer: you are plotting points singularly so you can't see them
short solution: pool them together in arrays then use the command to plot
arrays can be made using the line
[a,b,c,d,.....,z]
where letters are your values for bottom or top
I am not familiar with the sym data format but I think the issue might be just in the way you plot there. Also comment the keyboard otherwise it gets stuck in debugging.
see your code with scatter instead of plot (scatter can plot single points, plot can too but if you plot a single points one by one you will not see them, they have to be organized in arrays)
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;
% scatter(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};
scatter(b0(1,1), z(1,1))
% keyboard
hold on
scatter(t0(1,1), z(2,1))
scatter(b1(1,1), z(2,1))
scatter(t1(1,1), z(3,1))
scatter(b2(1,1), z(3,1))
scatter(t2(1,1), z(4,1))
scatter(b3(1,1), z(4,1))
scatter(t3(1,1), z(5,1))
scatter(b4(1,1), z(5,1))
scatter(t4(1,1), z(6,1))
scatter(b5(1,1), z(6,1))
scatter(t5(1,1), z(7,1))
scatter(b6(1,1), z(7,1))
scatter(t6(1,1), z(8,1))
scatter(b7(1,1), z(8,1))
scatter(t7(1,1), z(9,1))
scatter(b8(1,1), z(9,1))
scatter(t8(1,1), z(10,1))
scatter(b9(1,1), z(10,1))
scatter(t9(1,1), z(11,1))
scatter(b10(1,1), z(11,1))
scatter(t10(1,1), z(12,1))
scatter(b11(1,1), z(12,1))
scatter(t11(1,1), z(13,1))
scatter(b12(1,1), z(13,1))
scatter(t12(1,1), z(14,1))
scatter(b13(1,1), z(14,1))
scatter(t13(1,1), z(15,1))
scatter(b14(1,1), z(15,1))
scatter(t14(1,1), z(16,1))
scatter(b15(1,1), z(16,1))
scatter(t15(1,1), z(17,1))
hold off
% keyboard
0 comentarios
per isakson
el 17 de Abr. de 2021
I added 'd' to all your plot-commands
plot(b0(1,1), z(1,1), 'd' )
keyboard
hold on
plot(t0(1,1), z(2,1), 'd' )
plot(b1(1,1), z(2,1), 'd' )
plot(t1(1,1), z(3,1), 'd' )
plot(b2(1,1), z(3,1), 'd' )
etc.
Now the script outputs
3 comentarios
per isakson
el 17 de Abr. de 2021
Editada: per isakson
el 17 de Abr. de 2021
To draw a line you need a vector of at least two elements. That's true for both plot and stairs. I joined only a subset of your points.
x = [t0(1,1),b1(1,1),t1(1,1),b2(1,1),t2(1,1),b3(1,1),t3(1,1),b4(1,1),t4(1,1),b5(1,1),t5(1,1),b6(1,1),t6(1,1),b7(1,1),t7(1,1)];
y = [z(2,1),z(2,1),z(3,1),z(3,1),z(4,1),z(4,1),z(5,1),z(5,1),z(6,1),z(6,1),z(7,1),z(7,1),z(8,1),z(8,1),z(9,1)];
stairs( x, y, '-dk' )
outputs
and replace '-dk' by '-k' to remove the markers
Not all of the lines of unnamed.png are vertical or horizontal. plot( x, y, '-dk' ) produces
Ver también
Categorías
Más información sobre 2-D and 3-D Plots en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!