Deforming sheet into conical structure

I am trying to deform a flat sheet into a conical structure using MATLAB. (Pls see the image below)

L = length of sheet; R= bigger radius; r = smaller radius

Here's my approach:

1. First i create an upward ramp (as shown in the image above), with a displacement in z-direction given by:

u_1 = r/L * (pi-1) * x

2. Then curve the edges of rectangular sheet into the conical form, given by displacements:

3. Hence the net displacement of any point becomes u_1 + u_2.

I have written this as a matlab code:

for i=1:length(A)

A(i,2) = A(i,2) + (sqrt(2*A(i,3)*r - A(i,3)*A(i,3))-A(i,2))*(A(i,1)/L);

A(i,3) = A(i,3) + ((r/L)*(pi-1)*(A(i,1)) + (R-sqrt(R*R - A(i,2)*A(i,2))*(1-(A(i,1)/L))) + ((R-r)+ (2/pi)*A(i,2))*(A(i,1)/L));

A here is the matrix with initial coordinates of the flat sheet. Then using a for loop, i apply the u_1 + u_2 displacement, changing the y and z coordinates of the sheet.

However, after running the code, i get one side to be very nicely fitting the curvature of the cone, whereas on the other side, all the points coincide at one single point rather than forming a circular face. Please see the images below.

Original flat sheet

Curved conical sheet with the issue

Can someone please help me with this?

Answer 1

[William Rose]

@Priyanshu,

You want to map points from the x-y plane (flat sheet) to half of a conical frustrum.

The initial points are points (x1,y1) in sheet 1: x1=[0,L], y1=[-0.5,+0.5]. The final points are on half of a conical frustrum, with axis from (x2,y2,z2)=(0,0,R1) to (L,0,R1), and radius R1 at x2=0, radius R2 at x2=L. I have changed your notation somewhat, because I want to define r(x2) as a function.

The mapping from x1 to x2 is simple:

x2=x1 eq.0

The mapping from y1 to (y2,z2) is more complicated. Define a radius function:

r(x2)=R1+x2*(R2-R1)/L eq.1

Since x2=x1, we can write

r(x1)=R1+x1*(R2-R1)/L eq.2

Define angle theta as the angle measured about the cone axis, in a plane perpendicular to the axis, i.e. theta is the angle in the y-z plane. Define theta so that theta=0 when z2<0 and y2=0; theta=-pi/2 when y2<0 and z2=R; theta=pi/2 when y2>0 and z2=R.

It follows from the definition of r(x2) and theta that (r,theta) and (y2,z2) are related as follows:

y2=r(x1)*sin(theta) eq.3

z2=R1-r(x1)*cos(theta) eq.4

Now we choose to map y1 (domain=-0.5 to +0.5) to theta (range -pi/2 to +pi/2).

theta=pi*y1 eq.5

Substitute eq.2 and eq.5 into eq.3:

y2=(R1+x1*(R2-R1)/L)*sin(pi*y1) eq.6

Substitute eq.2 and eq.5 into eq.4:

z2=R1-(R1+x1*(R2-R1)/L)*cos(pi*y1) eq.7

Equations 0, 6, and 7 completely define the mapping from the sheet to the conical frustrum.

In Matlab:

L=1; R1=0.5; R2=0.2;
x1=L*(0:1/6:1); y1=-0.5:0.1:0.5;
[X1,Y1]=meshgrid(x1,y1);
X2=X1;
Y2=(R1+X1*(R2-R1)/L).*sin(pi*Y1);
Z2=R1-(R1+X1*(R2-R1)/L).*cos(pi*Y1);
% Plot results
figure;
colr=['k','r','g','b','c','m','y'];
subplot(211)
for i=1:7 
    plot3(X1(:,i),Y1(:,i),zeros(length(y1),1),'*','Color',colr(i)); hold on
end
axis equal; grid on; view(45,20); zlim([0 R1])
xlabel('X1'); ylabel('Y1'); zlabel('Z1')
subplot(212)
for i=1:7 
    plot3(X2(:,i),Y2(:,i),Z2(:,i),'*','Color',colr(i)); hold on
end
axis equal; grid on; view(45,20)
xlabel('X2'); ylabel('Y2'); zlabel('Z2')

It seems to work.

[Priyanshu]

@William Rose

Thank you for the help. The solution is easy to follow along and the explanation is great.

However, I am still facing a couple of issues:

1. The initial metal sheet is a 3-D strucuture, like a thin rectangular sheet with z:[0,1]. When I try including this in the solution, MATLAB is unable to update the z-coordinate.
2. Also for the mapping into conical part, the metal sheet should cover the lower half of the bigger face and then linearlty increase its curvature and finally cover the entire circumference of the smaller face. I beleive for this transformation, we'll need two theta values, one for the bigger face and one for the smaller face and then linearly increase the curvature of the flat sheet as a distance from either end.

I would be happy if you can let me know how to tackle these two issues.

[William Rose]

@Priyanshu, I have been having trouble accessing this site while travelling. I can answer your question. I hope to do so in the next day or two.

[William Rose]

@Priyanshu,

You want to map points from a flat sheet to part of a conical frustrum.

The initial points are points (x1,y1,z1) in sheet 1: x1=[0,L], y1=[-0.5,+0.5], z1=0. The final points are on half of a conical frustrum, with axis from (x2,y2,z2)=(0,0,R1) to (L,0,R1), and radius R1 at x2=0, radius R2 at x2=L. I have changed your notation somewhat, because I want to define r(x2) as a function.

The mapping from x1 to x2 is simple:

x2=x1 eq.0

The mapping from (y1,z1) to (y2,z2) is more complicated. Define a radius function:

r(x2)=R1+x2*(R2-R1)/L eq.1

Since x2=x1, we can write

r(x1)=R1+x1*(R2-R1)/L eq.2

Define theta as the angle measured about the cone axis, in a plane perpendicular to the axis, i.e. theta is the angle in the y-z plane. Define theta so that theta=0 when z2<0 and y2=0; theta=-pi/2 when y2<0 and z2=R; theta=pi/2 when y2>0 and z2=R.

It follows from the definition of r(x2) and theta that (r,theta) and (y2,z2) are related as follows:

y2=r(x1)*sin(theta) eq.3

z2=R1-r(x1)*cos(theta) eq.4

We want the flat sheet to wrap half-way around the cone at x1=0, and we want the flat sheet to wrap all the way around the cone at x1=L. In other words, thetaMax=pi/2 at x1=0, and thetaMax=pi at x`=L. Therefore we have

thetaMax(x1)=pi/2+x1*pi/(2*L) eq.5

and we map (x1,y1) to theta by

theta(x1,y1)=thetaMax(x1)*y1*2 eq.6

Substitute eq.2 and eq.6 into eq.3:

y2=(R1+x1*(R2-R1)/L)*sin(thetaMax*y1*2) eq.7

Substitute eq.5 into eq.7:

y2=(R1+x1*(R2-R1)/L)*sin((pi+pi*x1/L)*y1) eq.8

Substitute eq.2 and eq.6 into eq.4:

z2=R1-(R1+x1*(R2-R1)/L)*cos(thetaMax*y1*2) eq.9

Substitute eq.5 into eq.9:

z2=R1-(R1+x1*(R2-R1)/L)*cos((pi+pi*x1/L)*y1) eq.10

Equations 0, 8, and 10 completely define the mapping from the sheet to the conical frustrum.

In Matlab:

L=1; R1=0.5; R2=0.2;        % length, radii
nx=11; ny=13;               % number of x-values and y-values
x1=L*linspace(0,L,nx);      % x values
y1=0.5*linspace(-1,1,ny);   % y values
[X1,Y1]=meshgrid(x1,y1);
Z1=zeros(size(X1));
X2=X1;
Y2=(R1+X1*(R2-R1)/L).*sin((pi+pi*X1/L).*Y1);       % eq.8
Z2=R1-(R1+X1*(R2-R1)/L).*cos((pi+pi*X1/L).*Y1);    % eq.10
% Plot results
figure;
colr=hsv2rgb([linspace(0,1,nx)',ones(nx,2)]);
subplot(211)
for i=1:nx 
    plot3(X1(:,i),Y1(:,i),Z1(:,i),'*','Color',colr(i,:)); hold on
end
axis equal; grid on; view(45,20); zlim([0 R1])
xlabel('X1'); ylabel('Y1'); zlabel('Z1')
subplot(212)
for i=1:nx 
    plot3(X2(:,i),Y2(:,i),Z2(:,i),'*','Color',colr(i,:)); hold on
end
axis equal; grid on; view(45,20)
xlabel('X2'); ylabel('Y2'); zlabel('Z2')

OK

[Priyanshu]

@William Rose Thank you! This helps a lot. How about instead of initializing the metal flat sheet with z1=0, we also have an array of z1 values similar to the arrays of x and y values. And by using meshgrid function, we can make a 3D solid having x1,y1,z1. The diffferent values of z1 will give height/depth to the metal sheet.

[William Rose]

@Priyanshu,

In the most recent code, I did initialize Z as an array with the same dimensions as arrays X and Y.

Yes, you can use mesh() (not meshgrid()). In that case you would not use the for loop. I used a for loop above so that I could color each row separately.

L=1; R1=0.5; R2=0.2;        % length, radii
nx=11; ny=13;               % number of x-values and y-values
x1=L*linspace(0,L,nx);      % x values
y1=0.5*linspace(-1,1,ny);   % y values
[X1,Y1]=meshgrid(x1,y1);
Z1=zeros(size(X1));
X2=X1;
Y2=(R1+X1*(R2-R1)/L).*sin((pi+pi*X1/L).*Y1);       % eq.8
Z2=R1-(R1+X1*(R2-R1)/L).*cos((pi+pi*X1/L).*Y1);    % eq.10
% Plot results
figure;
colr=hsv2rgb([linspace(0,1,nx)',ones(nx,2)]);
subplot(211)
mesh(X1,Y1,Z1);
% for i=1:nx 
%     plot3(X1(:,i),Y1(:,i),Z1(:,i),'*','Color',colr(i,:)); hold on
% end
axis equal; grid on; view(45,20); zlim([0 R1])
xlabel('X1'); ylabel('Y1'); zlabel('Z1')
subplot(212)
mesh(X2,Y2,Z2)
% for i=1:nx 
%     plot3(X2(:,i),Y2(:,i),Z2(:,i),'*','Color',colr(i,:)); hold on
% end
axis equal; grid on; view(45,20)
xlabel('X2'); ylabel('Y2'); zlabel('Z2')

[Priyanshu]

@William Rose thank you!.

[William Rose]

@Priyanshu, you're welcome.

[Priyanshu]

@William Rose, In your solution, https://www.mathworks.com/matlabcentral/answers/2134556-deforming-sheet-into-conical-structure#comment_3207473

the z-height is varying sinosouidally, as shown in the image below.

I was trying to change your code to vary the z-heght linearly from X=0 to X=L, so that the maximum height occurs at X=L. I figured out that the sin and cos terms in Y and Z deformations play a role in this variation, however I am unable to acheive what I want.

Could you please guide me with this?

[William Rose]

@Priyanshu,

[Edit: Correct the last equation, at the end of this post.]

I did not expect that you would still be working on this.

In the most recent version of the sheet-wrapping algorithm, the initial sheet is specified by the set of points [X1, Y1, Z1], and the wrapped points are specified by the set [X2,Y2, Z2]. The wrapped points lie on the surface of a frustrum of a cone. The centerline of the cone is a line parallel to the x-axis, with Y2=0 and Z2=R1. This means the cross sections of the cone, perpendicular to its axis, are circles parallel to the Y2-Z2 plane. The circle radius when X2=0, and the circle radius when X2=L; i.e.

.

You requested that the wrapped surface extend half way around the cone at X2=0, and that the wrapped surface wrap all the way around the cone at X2=L. Therefore I made the total wrap angle, , be a linear function of X2: at X2=0, and at X2=L; i.e.

.

With these choices, the max height (max Z2 value) of the wrapped surface, for each value of X2, is a constant plus the product of a linearly increasing function and a sinusoid:

.

Now you want to be a linear function of x2. Is that right? Do you still want the points of the sheet to lie on the frustrum of a cone?

[William Rose]

@Priyanshu,

[Edit: change an equation for z2,max below.]

if you want the ponts of the sheet to lie on the frustrum of a cone, then you cannot change the function , because it represents part of the definition of a cone. Therefore you must think about changing the function , in such a way that

increases linearly with x2. It follows, from the equations in my previous comment, that since R2<R1, is a linearly DEcreasing function of x2.

We want to increase linearly, from at X2=0 to at X2=L. Therefore we want

.

Therefore

That is my proposed new function for the total wrap-around angle . I suggest you modify the code I provided earlier to use this result. See if it gives a height that increases linearly with x2, as you requested.

[Priyanshu]

Hii @William Rose, I modified the Z2 to:

X2 = X1
Y2 = (radius1 + (radius2 - radius1) .* X / length) .* sin((pi + pi * X / length) .* Y);
Z2 = radius1 + (radius1 + (radius2 - radius1) .* X / length) .* sin (((2 * asin((radius2 .* X / length) / (radius1 + (radius2 - radius1) .* X / length)) + pi) - pi) / 2);

including the equations you gave above. However, I am unable to see the deformation. Here's the output figure:

[William Rose]

@Priyanshu,

Please post a complete script, and run it in the Matlab Answers window, so others can see what results or errors you get. Three isolated lines is not enough to understand the context. Thank you.

[Priyanshu]

@William Rose, apologies for the poorly worded comment. Here's the fullscript.

clc; clearvars; close all;
% Define constants
length = 1;
radius1 = length / (2);
radius2 = 1 * length / (2 * pi);

% Define grid dimensions
numX = 10; 
numY = 10;

% Generate grid points
xValues = linspace(0, length, numX) * length;
yValues = linspace(-1, 1, numY) * 0.5; 

% Create mesh grid
[X, Y] = meshgrid(xValues, yValues);
Z = zeros(size(X)); % Initial Z values are zero

% Creating matrix of intial points
X_init = X'; X_init = reshape(X_init, [numel(X),1]);
Y_init = Y'; Y_init = reshape(Y_init, [numel(Y),1]);
Z_init = Z'; Z_init = reshape(Z_init, [numel(Z),1]);
Init_points = [X_init, Y_init, Z_init];

% Writing the intial points into a text file for further calculations
writematrix(Init_points, 'init_test.txt','Delimiter', 'tab');

% Adjust Y and Z for deformation
Y_deformed = (radius1 + (radius2 - radius1) .* X / length) .* sin((pi + pi * X / length) .* Y);
Z_deformed = radius1 + (radius1 + (radius2 - radius1) .* X / length) .* sin (((2 * asin((radius2 .* X / length) / (radius1 + (radius2 - radius1) .* X / length)) + pi) - pi) / 2);

Warning: Matrix is singular to working precision.

% Calulating displacement vectos and final displacement matrix
X_disp = X - X; X_disp = X_disp'; X_disp = reshape(X_disp, [numel(X),1]);
Y_disp = Y_deformed - Y; Y_disp = Y_disp'; Y_disp = reshape(Y_disp, [numel(Y),1]);
Z_disp = Z_deformed - Z; Z_disp = Z_disp'; Z_disp = reshape(Z_disp, [numel(Z),1]);
Disp_points = [X_disp, Y_disp, Z_disp];


% Writing the displaced points into a text file for further calculations
writematrix(Disp_points, 'def_test.txt','Delimiter', 'tab');

% Calculating net displacement magnitudes
Displacement = sqrt((X - X).^2 + (Y_deformed - Y).^2 + (Z_deformed - Z).^2);

% Initialize figure
figure;

% Plot initial shape
subplot(2, 2, 1)
mesh(X, Y, Z);
axis equal; grid on; view(45, 20); zlim([0 radius1])
title('Initial Shape')
xlabel('X'); ylabel('Y'); zlabel('Z');

% Plot deformed shape colored by displacement
subplot(2, 2, 2)
mesh(X, Y_deformed, Z_deformed, Displacement, 'FaceAlpha','0.5', 'FaceColor','interp');% Add displacement as color data
colorbar % Show color bar to indicate scale of displacement
colormap jet % Choose a colormap that spans a wide range of visible colors
axis equal; grid on; view(45, 20)
xlabel('X'); ylabel('Y'); zlabel('Z');
title('Deformed Shape Colored by Displacement Magnitude')

[William Rose]

@Priyanshu, thank you. I will look at it.

[William Rose]

@Priyanshu,

From above, the full range for theta, at each value of x2, is:

% Define constants
L=1; R1=0.5; R2=0.2;        % length, radii
nx=11; ny=13;               % number of x-values and y-values
x1=L*linspace(0,L,nx);      % x values
y1=0.5*linspace(-1,1,ny);   % y values
% Make sheet 1
[X1,Y1]=meshgrid(x1,y1);
Z1=zeros(size(X1));
% Make sheet 2
X2=X1;
thetaFR=2*asin((x1*R2/L)./(R1+x1*(R2-R1)/L))+pi; %theta(full range)
theta=y1'*thetaFR;
Y2=(R1+X1*(R2-R1)/L).*sin(theta);
Z2=R1-(R1+X1*(R2-R1)/L).*cos(theta);

% Plot results
figure;
subplot(121)
mesh(X1,Y1,Z1);
axis equal; grid on; view(45,20); zlim([0 R1+2*R2])
xlabel('X1'); ylabel('Y1'); zlabel('Z1')
subplot(122)
mesh(X2,Y2,Z2)
axis equal; grid on; view(45,20); zlim([0 R1+2*R2])
xlabel('X2'); ylabel('Y2'); zlabel('Z2'); colormap('hsv')

If you click and drag the plot of the frustrum of a cone, you will see that the max z-height increases linearly with X2. That is what you requested. See screenshot below.

[Priyanshu]

Thank you @William Rose!

[William Rose]

@Priyanshu, you're welcome.