PDE toolbox: unexpected result of static structural model for stretchy sheet

I'm trying to model spatial deformation of a triangular stretchy sheet of material (think stretchy fabric or rubber). The top corner is fixed in place, while the bottom two corners at pulled with a specified force as shown below. Imagine holding one bottom vertex in each hand and pulling to make the sheet taut. I'd like to solve for the final shape of the sheet.

The results seems completely wrong, which means I am doing something completely wrong. I expect to see displacement of the x coordinate (result.Displacement.ux) change from negative to positive as x increases. Instead, x displacement mostly increases with y. If I decrease the mesh size the result changes apparently randomly. Any hint what I am doing wrong? Code below.

Code:

%% Generate x,y coordinates of edges of a triangle
a.XY1 = [0 0; 1.5 0; .75 3];            % (m)

% Create geometry in pde toolbox form
gd = [2; size(a.XY1,1); a.XY1(:,1); a.XY1(:,2)];
sf = 'P1';
dl = decsg(gd,sf,double(sf(:)));

%% Create model and apply geometry
model = createpde('Structural','static-planestrain');
geometryFromEdges(model, dl);

%% Add structural properties 
structuralProperties(model,'YoungsModulus',1e9,'PoissonsRatio',.4); % (Pa) (unitless)

% Add boundary conditions
structuralBC(model, 'Vertex',  3,'Constraint','fixed');

% Add load
f = 1e3; % (N)
structuralBoundaryLoad(model,'Vertex',1,'Force',[-1 -1]*f); 
structuralBoundaryLoad(model,'Vertex',2,'Force',[1 -1]*f);  

% Add mesh
generateMesh(model); 

% Solve
result = solve(model)
fprintf('Max x: %.2g, max y: %.2g\n',max(result.Displacement.ux),max(result.Displacement.uy))

%% Plot stuff
figure;
subplot(221)
pdegplot(dl,'VertexLabels','on') % Plot PDE geometry
title('Geometry'); axis off

subplot(222); 
pdeplot(model) 
title('Mesh')

subplot(223)
pdeplot(model,'XYData',result.Displacement.ux)
title('x-displacement')
axis equal

subplot(224)
pdeplot(model,'XYData',result.Displacement.uy)
title('y-displacement')
colormap('jet')
axis equal

Answer 1

[Pratyush Swain]

Hi Myrtle42,

One key observation which can be noted here is the forces are applied with the assumption that force vectors are [-1,-1] and [1,-1] scaled by the force magnitude f=1e3.This implies that the forces are being applied at a 45-degree angle downwards from the horizontal, which would be correct if the triangle was a right isosceles triangle with the right angle at the top vertex. However, in the defined triangle, the vertices are defined as [0 0; 1.5 0;.75 3], which does not form a right isosceles triangle.

Please refer to this revised code which calculates the actual direction of the forces based on the geometry of the given triangle.This is done by finding the unit vectors along the edges from the fixed vertex to the vertices where the forces are applied.

%% Generate x,y coordinates of edges of a triangle
a.XY1 = [0 0; 1.5 0; .75 3];            % (m)

%% Create geometry in pde toolbox form
gd = [2; size(a.XY1,1); a.XY1(:,1); a.XY1(:,2)];
sf = 'P1';
dl = decsg(gd,sf,double(sf(:)));

%% Create model and apply geometry
model = createpde('Structural','static-planestrain');
geometryFromEdges(model, dl);

%% Add structural properties 
structuralProperties(model,'YoungsModulus',1e9,'PoissonsRatio',.4);

%% Add boundary conditions
structuralBC(model, 'Vertex',  3,'Constraint','fixed');

%% Coordinates of vertices
v1 = a.XY1(1,:);
v2 = a.XY1(2,:);
v3 = a.XY1(3,:);

%% Calculate unit vectors along the edges from v3 to v1 and v3 to v2
%% This ensures that the forces are applied precisely along the direction of the edges
%% Important for an accurate simulation of the physical scenario described.
edge1 = v1 - v3;
edge2 = v2 - v3;
unitVector1 = edge1 / norm(edge1);
unitVector2 = edge2 / norm(edge2);


%% Force magnitude
f = 1e3; % (N)

%% Apply forces in the direction of the edges
force1 = unitVector1 * f;
force2 = unitVector2 * f;

%% Apply the calculated forces to the model
structuralBoundaryLoad(model,'Vertex',1,'Force',force1); 
structuralBoundaryLoad(model,'Vertex',2,'Force',force2);

%% Add mesh with specified maximum element size
generateMesh(model, 'Hmax', 0.1); 

%% Solve
result = solve(model);
fprintf('Max x: %.2g, max y: %.2g\n',max(result.Displacement.ux),max(result.Displacement.uy));

Max x: 1.2e-05, max y: 0

%% Plot stuff
figure;
subplot(221)
pdegplot(dl,'VertexLabels','on') % Plot PDE geometry
title('Geometry'); axis equal; axis off;
subplot(222); 
pdeplot(model) 
title('Mesh')
subplot(223)
pdeplot(model,'XYData',result.Displacement.ux, 'Deformation',result.Displacement,'DeformationScaleFactor',1)
title('x-displacement')
axis equal
subplot(224)
pdeplot(model,'XYData',result.Displacement.uy, 'Deformation',result.Displacement,'DeformationScaleFactor',1)
title('y-displacement')
colormap('jet')
axis equal

For more information on pde toolbox please refer to https://www.mathworks.com/help/pde/

Hope this helps.