I face a problem with the dimensions of Z as it should be a 2x2 matrix.
62 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Mohd Babar
el 30 de Jul. de 2024 a las 5:25
Comentada: Walter Roberson
el 30 de Jul. de 2024 a las 8:29
%------Creating a strain matrix-----------%
E = zeros(2*(size(K1,1)-1),2*(size(K1,1)-1));
for j = 1: (size(K1,1)-1)
for i = 1: (size(K1,1)-1)
v = [zpk1(j,i);zpk2(j,i);zpk1(j,i+1);zpk2(j,i+1);zpk1(j+1,i+1);zpk2(j+1,i+1);zpk1(j+1,i);zpk2(j+1,i)];
e1= strain1(v);
e2= strain2(v);
e3= strain3(v);
e4= strain4(v);
E((2*j-1),(2*i-1)) = e1(1);
E((2*j-1),2*i) = e2(1);
E(2*j,(2*i-1)) = e4(1);
E(2*j,2*i) = e3(1);
end
end
su = round(4*d/g)+1;
r = E((((su-1)- round(d/g)):-1:1),(su-1));
p = E((((su-1)- round(d/g)):-1:1),(su));
st = (p+r)./(1.5*(10)^(-4));
m = g*(size(r)-1);
H = 0:g:(m);
h = H./10;
writematrix(st)
type 'st.txt'
writematrix(h)
type 'h.txt'
plot(h,st)
xlim([0 1.5])
xlabel("X2/d",'fontsize',14)
ylabel("Normalised Strain along vertical path", 'fontsize',13)
legend("d/w = 0.5")
saveas(gcf,'plot.png')
contourf(st)
colorbar
saveas(gcf,'contour.png')
The final plot of the above code should resemble something like this:
I have double chekced, the formula is correct and the value of st at a given r and p match with the experimental results. However, I face a problem with the dimensions of Z as it should be a 2x2 matrix.
Thank you in advance.
3 comentarios
Mohd Babar
el 30 de Jul. de 2024 a las 6:24
Editada: Walter Roberson
el 30 de Jul. de 2024 a las 6:27
Respuestas (1)
Walter Roberson
el 30 de Jul. de 2024 a las 6:36
d = 10;
g = 0.085;
d and g are scalars
su = round(4 * d / g) + 1;
su is formed from scalar d and g, so su is scalar
r = E((((su - 1) - round(d / g)):-1:1), (su - 1));
p = E((((su - 1) - round(d / g)):-1:1), su);
su-1 and su are scalars, so you are extracting exactly one column from E into r and p, so r and p are vectors (not 2D)
st = (p + r) / (1.5 * 10^(-4));
vector plus vector (in the same direction) is vector, so st is a vector.
contourf(st)
You are trying to contourf() a vector.
5 comentarios
Walter Roberson
el 30 de Jul. de 2024 a las 8:21
Editada: Walter Roberson
el 30 de Jul. de 2024 a las 8:28
Using meshgrid like that is not going to work. You have
st = (p + r) / (1.5 * 10^(-4));
and making p and r into grids cannot possibly give you the kind of curved output that you hope for.
clc;
clear all;
%-----Geometry, Number of simulations------%
d = 10;
N = 1;
a11 = 250 * (d / 10);
a12 = 100 * (d / 10);
%-----Mesh grid boundary definition--------%
xms = a11 / 2 - 2 * d;
xmf = a11 / 2 + 2 * d;
yms = a12 / 2 - 2 * d;
ymf = a12 / 2 + 2 * d;
%----------------Grid size----------------%
g = 0.085;
xq1 = xms:g:xmf;
xq = xq1';
yq1 = yms:g:ymf;
yq = yq1';
%-----Mesh grid around the hole----------%
[K1, K2] = meshgrid(xq, yq);
a = size(xq1, 2);
b = size(yq1, 2);
for k = 1:b
for j = 1:a
e = ((k * g) - 2 * d)^2 + ((j * g) - 2 * d)^2 - (0.5 * d)^2;
if e < 0
K1(j, k) = NaN;
K2(j, k) = NaN;
end
end
end
zz1 = 0;
zz2 = 0;
%----Importing files from Abaqus---------%
for i = 1:N
fx = ['aax' num2str(i) '.txt'];
fy = ['aay' num2str(i) '.txt'];
fu = ['aau1' num2str(i) '.txt'];
fv = ['aau2' num2str(i) '.txt'];
x = textfile(fx);
y = textfile(fy);
z1 = textfile(fu);
z2 = textfile(fv);
U1 = griddata(x, y, z1, K1, K2);
U2 = griddata(x, y, z2, K1, K2);
zz1 = U1 + zz1;
zz2 = U2 + zz2;
end
%--------Averaging the displacements-------%
zpk1 = zz1 / N;
zpk2 = zz2 / N;
xlswrite('zpk1.xlsx', zpk1)
xlswrite('zpk2.xlsx', zpk2)
%------Creating a strain matrix-----------%
E = zeros(2 * (size(K1, 1) - 1), 2 * (size(K1, 1) - 1));
for j = 1:(size(K1, 1) - 1)
for i = 1:(size(K1, 1) - 1)
v = [zpk1(j, i); zpk2(j, i); zpk1(j, i + 1); zpk2(j, i + 1); zpk1(j + 1, i + 1); zpk2(j + 1, i + 1); zpk1(j + 1, i); zpk2(j + 1, i)];
e1 = strain1(v);
e2 = strain2(v);
e3 = strain3(v);
e4 = strain4(v);
E((2 * j - 1), (2 * i - 1)) = e1(1);
E((2 * j - 1), 2 * i) = e2(1);
E(2 * j, (2 * i - 1)) = e4(1);
E(2 * j, 2 * i) = e3(1);
end
end
su = round(4 * d / g) + 1;
r = E((((su - 1) - round(d / g)):-1:1), (su - 1));
p = E((((su - 1) - round(d / g)):-1:1), su);
[r, p] = meshgrid(r,p);
st = (p + r) / (1.5 * 10^(-4));
m = g * (size(r, 1) - 1);
H = 0:g:m;
h = H / 10;
writematrix(st, 'st.txt')
writematrix(h, 'h.txt')
figure;
plot(h, st)
xlim([0 1.5])
xlabel("X2/d", 'fontsize', 14)
ylabel("Normalized Strain along vertical path", 'fontsize', 13)
legend("d/w = 0.5")
saveas(gcf, 'plot.png')
figure;
contourf(st)
colorbar
saveas(gcf, 'contour.png')
function E1 = strain1(v)
B1 = [1 0 0 0; 0 0 0 1; 0 1 1 0];
x1 = -1 / sqrt(3); eta1 = -1 / sqrt(3);
B21 = B2();
B31 = B3(x1, eta1);
B12 = B1 * B21;
B = B12 * B31;
E1 = B * v;
end
function E1 = strain2(v)
B1 = [1 0 0 0; 0 0 0 1; 0 1 1 0];
x1 = 1 / sqrt(3); eta1 = -1 / sqrt(3);
B21 = B2();
B31 = B3(x1, eta1);
B12 = B1 * B21;
B = B12 * B31;
E1 = B * v;
end
function E1 = strain3(v)
B1 = [1 0 0 0; 0 0 0 1; 0 1 1 0];
x1 = 1 / sqrt(3); eta1 = 1 / sqrt(3);
B21 = B2();
B31 = B3(x1, eta1);
B12 = B1 * B21;
B = B12 * B31;
E1 = B * v;
end
function E1 = strain4(v)
B1 = [1 0 0 0; 0 0 0 1; 0 1 1 0];
x1 = -1 / sqrt(3); eta1 = 1 / sqrt(3);
B21 = B2();
B31 = B3(x1, eta1);
B12 = B1 * B21;
B = B12 * B31;
E1 = B * v;
end
function B21 = B2()
B21 = [0.05 / 0.0025 0 0 0; 0 0.05 / 0.0025 0 0; 0 0 0.05 / 0.0025 0; 0 0 0 0.05 / 0.0025];
end
function B31 = B3(x1, eta1)
B31 = [...
-(1 - eta1) * 0.25 0 (1 - eta1) * 0.25 0 (1 + eta1) * 0.25 0 -(1 + eta1) * 0.25 0; ...
-(1 - x1) * 0.25 0 -(1 + x1) * 0.25 0 (1 + x1) * 0.25 0 (1 - x1) * 0.25 0; ...
0 -(1 - eta1) * 0.25 0 (1 - eta1) * 0.25 0 (1 + eta1) * 0.25 0 -(1 + eta1) * 0.25; ...
0 -(1 - x1) * 0.25 0 -(1 + x1) * 0.25 0 (1 + x1) * 0.25 0 (1 - x1) * 0.25];
end
function value = textfile(filename)
fileID = fopen(filename, 'r');
if fileID == -1
error('Failed to open file: %s', filename);
end
value = fscanf(fileID, '%f');
fclose(fileID);
end
Walter Roberson
el 30 de Jul. de 2024 a las 8:29
Ver también
Categorías
Más información sobre Bessel functions 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!