Note
Go to the end to download the full example code.
Elas7#
Control lever for a molding machine used to blow plastic bottles.
References#
This example comes from:
A French research article: MMC et RDM comparées sur un cas de dimensionnement de levier
The book: Mécanique des systèmes et des milieux déformables (3rd edition)


==================== Mesh ====================
Element type: PRISM6
Ne = 6108, Nn = 4315
==================== Model ====================
Isotropic:
E = 2.10e+05, v = 0.25
solver : scipy
============= Boundary Conditions =============
Unspecified.
=================== Results ===================
W def = 615.59
Svm max = 43.45
Evm max = 0.03 %
Ux max = 2.07e-02
Ux min = -2.05e-02
Uy max = 1.98e-03
Uy min = -1.13e-01
Uz max = 4.60e-04
Uz min = -4.60e-04
=================== TicTac ===================
Mesh : 323.721 ms
Boundary Conditions : 167.608 µs
Matrix : 688.181 ms
Solver : 3.332 s
PostProcessing : 17.106 ms
Generate paraview 1/5 (20.00 %) 278.38 ms
Generate paraview 2/5 (40.00 %) 208.27 ms
Generate paraview 3/5 (60.00 %) 134.76 ms
Generate paraview 4/5 (80.00 %) 68.33 ms
Generate paraview 5/5 (100.00 %) 0.00 µs
20 import numpy as np
21
22 from EasyFEA import Display, Folder, Models, ElemType, Simulations, PyVista, Paraview
23 from EasyFEA.Geoms import Points, Point, Circle
24
25 if __name__ == "__main__":
26 Display.Clear()
27
28 # ----------------------------------------------
29 # Configuration
30 # ----------------------------------------------
31
32 # outputs
33 folder = Folder.Results_Dir()
34 makeMovie = False
35 makeParaview = True
36
37 # geom
38 dim = 3
39 thickness = 25
40
41 # sinusoidal load
42 F = 13000 # N
43 area = np.pi * 25 / 2 * thickness
44 coefs = np.linspace(0, 1, 5)
45 loads = F * np.sin(np.pi / 2 * coefs)
46
47 # ----------------------------------------------
48 # Mesh
49 # ----------------------------------------------
50 meshSize = 42 / 5
51
52 # get the contour
53 pt1 = Point(0, 80)
54 pt2 = Point(-245, 21)
55 pt3 = Point(-245, -21)
56 pt4 = Point(0, -80)
57 pt5 = Point(114, -80)
58 pt6 = Point(114 + 13, -80 + 13)
59 pt7 = Point(114 + 13, -67 / 2)
60 pt8 = Point(114 + 13 + 25, -67 / 2)
61 pt9 = Point(114 + 13 + 25, -34 / 2)
62 pt10 = Point(114 + 13, -34 / 2, r=10)
63 pt11 = Point(114 + 13, 34 / 2, r=10)
64 pt12 = Point(114 + 13 + 25, 34 / 2)
65 pt13 = Point(114 + 13 + 25, 67 / 2)
66 pt14 = Point(114 + 13, 67 / 2)
67 pt15 = Point(114 + 13, 80 - 13)
68 pt16 = Point(114, 80)
69
70 # fmt: off
71 contour = Points([
72 pt1, pt2, pt3, pt4, pt5, pt6, pt7, pt8, pt9, pt10,
73 pt11, pt12, pt13, pt14, pt15, pt16,
74 ],
75 meshSize,
76 )
77 # fmt: on
78
79 # get circles
80 circle1 = Circle(Point(-220), 25, meshSize, isHollow=True)
81 circle2 = Circle(Point(), 25, meshSize, isHollow=True)
82 circle3 = Circle(Point(114 + 13 - 25, 80 - 25), 25, meshSize, isHollow=True)
83 circle4 = Circle(Point(114 + 13 - 25, -80 + 25), 25, meshSize, isHollow=True)
84 circle5 = Circle(Point(67), 84, meshSize, isHollow=True)
85
86 inclusions = [circle1, circle2, circle3, circle4, circle5]
87
88 # get the mesh
89 if dim == 2:
90 mesh = contour.Mesh_2D(inclusions, ElemType.TRI3)
91 else:
92 mesh = contour.Mesh_Extrude(inclusions, [0, 0, thickness], [4], ElemType.PRISM6)
93
94 # get loading nodes
95 nodesLoad = mesh.Nodes_Cylinder(circle1)
96 nodesLoad = nodesLoad[mesh.coord[nodesLoad, 1] <= 0]
97
98 # get fixed nodes
99 nodesCircle2 = mesh.Nodes_Cylinder(circle2)
100 nodesCircle3 = mesh.Nodes_Cylinder(circle3)
101 nodesCircle4 = mesh.Nodes_Cylinder(circle4)
102 fixedNodes = np.concat((nodesCircle2, nodesCircle3, nodesCircle4))
103
104 # PyVista.Plot_Mesh(mesh).show()
105
106 # ----------------------------------------------
107 # Simulation
108 # ----------------------------------------------
109
110 material = Models.Elastic.Isotropic(
111 dim, E=210000, v=0.25, planeStress=True, thickness=thickness
112 )
113
114 simu = Simulations.Elastic(mesh, material)
115
116 fixed = [0] * dim
117 unknowns = simu.Get_unknowns()
118
119 # loop over the load
120 for load in loads:
121 # Apply boundary conditions
122 simu.Bc_Init()
123 simu.add_dirichlet(fixedNodes, fixed, unknowns)
124 simu.add_surfLoad(nodesLoad, [-load / area], ["y"])
125
126 # Solve the simulation
127 simu.Solve()
128
129 # Save the iteration results
130 simu.Save_Iter()
131
132 # ----------------------------------------------
133 # Results
134 # ----------------------------------------------
135
136 print(simu)
137
138 PyVista.Plot_BoundaryConditions(simu).show()
139
140 PyVista.Plot(simu, "Svm", deformFactor=200, plotMesh=True, nodeValues=False).show()
141
142 if makeMovie:
143 PyVista.Movie_simu(
144 simu,
145 "Svm",
146 folder,
147 "Svm.gif",
148 deformFactor=200,
149 nodeValues=False,
150 plotMesh=True,
151 )
152
153 if makeParaview:
154 Paraview.Save_simu(
155 simu, folder, elementFields=["Svm", "Stress", "Strain", "Wdef_e"]
156 )
Total running time of the script: (0 minutes 6.134 seconds)