Skip to content

Deflection Analysis of Bracket¤

This example shows how to analyze a 3-D mechanical part under an applied load using the finite element analysis model and determine the maximal deflection.

Matlab examples see Partial Differential Equation Toolbox.

Partial Differential Equation Toolbox Solving¤

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
clc; clear;


E  = 200e9;
nu = 0.3;
rho = 0.0;


model = femodel(AnalysisType="structuralStatic", ...
                Geometry="BracketWithHole.stl");


model.MaterialProperties = materialProperties( ...
    YoungsModulus = E, ...
    PoissonsRatio = nu);


model.FaceBC(4) = faceBC(Constraint="fixed");
model.FaceLoad(8) = faceLoad(SurfaceTraction=[0;0;-1e4]);


figure
pdegplot(model,FaceLabels="on");
view(30,30);
title("Bracket with Face Labels")
figure_0.png
1
2
3
4
5
6
pdeModel = generateMesh(model, GeometricOrder="linear");


figure
pdemesh(pdeModel);
title("Mesh with Linear Tetrahedral Elements")
figure_1.png
1
2
3
result = solve(pdeModel);
minUzPDE = min(result.Displacement.uz);
fprintf("Maximal deflection in the z-direction is %g meters.",minUzPDE)
Output
Maximal deflection in the z-direction is -1.44604e-05 meters. 
1
2
3
4
figure
pdeplot3D(result.Mesh,ColorMapData=result.Displacement.uz);
title("z-displacement");
colormap("jet")
figure_2.png
1
[az, el] = view;

1
2
3
4
figure
pdeplot3D(result.Mesh,ColorMapData=result.VonMisesStress)
title("von Mises stress")
colormap("jet")
figure_3.png
1
2
maxVonMisesStressPDE = max(result.VonMisesStress);
fprintf("Maximal VonMisesStress is %g Pa.", maxVonMisesStressPDE)
Output
Maximal VonMisesStress is 1.03705e+06 Pa.

OpenSeesMatlab Solving¤

Extract PDE mesh data¤

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
mesh = pdeModel.Mesh;


% PDE toolbox:
%   mesh.Nodes    : 3 x nNode
%   mesh.Elements : nNodePerElem x nElem
nodes = mesh.Nodes.';        % nNode x 3
elems = mesh.Elements.';     % nElem x nen


nNode = size(nodes,1);
nElem = size(elems,1);
nen   = size(elems,2);


if nen ~= 4
    error(['This script currently expects linear 4-node tetrahedra. ', ...
           'Please use generateMesh(..., GeometricOrder="linear").']);
end


%% ------------------------------------------------------------------------
% Find face nodes for fixed boundary and loaded boundary
%% ------------------------------------------------------------------------
fixedFaceID  = 4;
loadedFaceID = 8;
traction     = [0; 0; -1e4];   % N/m^2


fixedNodeIDs = findNodes(mesh, "region", "Face", fixedFaceID);
loadNodeIDs  = findNodes(mesh, "region", "Face", loadedFaceID);


fixedNodeIDs = unique(fixedNodeIDs(:));
loadNodeIDs  = unique(loadNodeIDs(:));


%% ------------------------------------------------------------------------
% Build boundary triangles on the loaded face and convert traction
%    to equivalent nodal loads
%% ------------------------------------------------------------------------
% For a tetra mesh, each tetra has 4 triangular faces.
% Boundary faces are those triangle faces that appear only once.
allBoundaryTris = getBoundaryTrianglesFromTets(elems);


% Keep only boundary triangles whose 3 nodes are all on Face 8
isOnLoadFace = all(ismember(allBoundaryTris, loadNodeIDs.'), 2);
loadFaceTris = allBoundaryTris(isOnLoadFace, :);


if isempty(loadFaceTris)
    warning('No boundary triangles were found on loaded Face %d.', loadedFaceID);
end


% Assemble nodal force vector from surface traction
nodalLoads = zeros(nNode, 3);


for i = 1:size(loadFaceTris,1)
    tri = loadFaceTris(i,:);
    xyz = nodes(tri, :);              % 3 x 3
    A   = triangleArea3D(xyz);


    % Constant traction on linear triangle:
    % equivalent nodal load = traction * A / 3 for each of 3 nodes
    fe = (A / 3) * traction(:);       % 3 x 1 for one node
    nodalLoads(tri, :) = nodalLoads(tri, :) + repmat(fe.', 3, 1);
end

Build OpenSees model¤

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
opsMAT = OpenSeesMatlab();
ops = opsMAT.opensees;


ops.wipe();
ops.model('basic', '-ndm', 3, '-ndf', 3);


% Nodes
for i = 1:nNode
    ops.node(i, nodes(i,1), nodes(i,2), nodes(i,3));
end


% Material
matTag = 1;
E  = 200e9;
nu = 0.3;
rho = 0.0;


% 3D isotropic elastic material
ops.nDMaterial('ElasticIsotropic', matTag, E, nu, rho);


% Elements: 4-node tetrahedra
for e = 1:nElem
    n = elems(e,:);
    ops.element('FourNodeTetrahedron', e, n(1), n(2), n(3), n(4), matTag);
end


% Apply fixed boundary from PDE Face 4
for k = 1:numel(fixedNodeIDs)
    nd = fixedNodeIDs(k);
    ops.fix(nd, 1, 1, 1);
end


% Apply equivalent nodal loads from PDE Face 8 traction
ops.timeSeries('Linear', 1);
ops.pattern('Plain', 1, 1);


tol = 0;
for i = 1:nNode
    fx = nodalLoads(i,1);
    fy = nodalLoads(i,2);
    fz = nodalLoads(i,3);


    if abs(fx) > tol || abs(fy) > tol || abs(fz) > tol
        ops.load(i, fx, fy, fz);
    end
end

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
% Static analysis
Nsteps = 1;
ops.constraints('Plain');
ops.numberer('RCM');
ops.system('BandGeneral');
ops.test('NormDispIncr', 1e-8, 50);
ops.algorithm('Linear');
ops.integrator('LoadControl', 1.0/Nsteps);
ops.analysis('Static');


ODB = opsMAT.post.createODB("myODB");  % create ODB
ok = ops.analyze(Nsteps);

Visualization by OpenSeesMatlab¤

1
2
3
opts = opsMAT.vis.defaultPlotModelOptions;
opts.loads.showNodal=true;
opsMAT.vis.plotModel(opts=opts);
Output
[OpenSeesMatlab] Model summary Nodes: 1336 Solid elements: 3883 
1
view(az, el);
figure_4.png
1
2
3
4
5
6
7
nodeResp = opsMAT.post.getNodalResponse("myODB");
uz = nodeResp.disp.uz;
minUz = min(uz(:));
fprintf("Maximal deflection in the z-direction:\n" + ...
        "  OpenSeesMatlab: %g meters\n" + ...
        "  MATLAB PDE    : %g meters\n", ...
        minUz, minUzPDE);
Output
Maximal deflection in the z-direction: OpenSeesMatlab: -1.44604e-05 meters MATLAB PDE : -1.44604e-05 meters 
1
2
3
4
5
6
7
8
9
opts = opsMAT.vis.defaultPlotNodalResponseOptions;
opts.fixed.show = false;
opts.surf.showEdges = false;


opsMAT.vis.plotNodalResponse(nodeResp, respType="disp", stepIdx="absMax", respComponent="UZ", opts=opts);
colormap("jet")
axis off
view(az, el);
figure_5.png
1
2
3
4
5
6
7
solidResp = opsMAT.post.getElementResponse("myODB", eleType="Solid");
vm = solidResp.StressMeasureAtNode.vonMises;
maxVonMisesStress = max(vm(:));
fprintf("Maximal Von Mises stress:\n" + ...
        "  OpenSeesMatlab: %g Pa\n" + ...
        "  MATLAB PDE    : %g Pa\n", ...
        maxVonMisesStress, maxVonMisesStressPDE);
Output
Maximal Von Mises stress: OpenSeesMatlab: 1.03705e+06 Pa MATLAB PDE : 1.03705e+06 Pa 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
opts = opsMAT.vis.defaultPlotContinuumResponseOptions;
opts.fixed.show = false;
opts.surf.showEdges = false;


opsMAT.vis.plotContinuumResponse(solidResp, stepIdx="absMax",...
    respType="StressMeasureAtNode", respComponent="vonMises", opts=opts);
axis off
colormap("jet")
view(az, el);
figure_6.png
1
% opsMAT.post.writeResponsePVD("myODB");

Local functions¤

1
2
3
%% ------------------------------------------------------------------------
% Local functions
%% ------------------------------------------------------------------------

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
function btris = getBoundaryTrianglesFromTets(tets)
%GETBOUNDARYTRIANGLESFROMTETS Return unique boundary triangle faces
% tets: nElem x 4


    f1 = tets(:, [1 2 3]);
    f2 = tets(:, [1 2 4]);
    f3 = tets(:, [1 3 4]);
    f4 = tets(:, [2 3 4]);


    allFaces = [f1; f2; f3; f4];
    allFacesSorted = sort(allFaces, 2);


    [uFaces, ~, ic] = unique(allFacesSorted, 'rows');
    counts = accumarray(ic, 1);


    isBoundary = counts(ic) == 1;
    boundaryFacesSorted = allFacesSorted(isBoundary, :);


    % Keep orientation-free representation
    btris = boundaryFacesSorted;
end


function A = triangleArea3D(xyz)
%TRIANGLEAREA3D Area of a 3D triangle
% xyz: 3 x 3, each row is one node coordinate


    v1 = xyz(2,:) - xyz(1,:);
    v2 = xyz(3,:) - xyz(1,:);
    A = 0.5 * norm(cross(v1, v2));
end