Note
Go to the end to download the full example code.
Elas1#
A cantilever beam undergoing bending deformation.
==================== Mesh ====================
Element type: QUAD9
Ne = 81, Nn = 385
==================== Model ====================
ElasIsot:
E = 2.10e+05, v = 0.3
planeStress = True
thickness = 1.30e+01
solver : scipy
============= Boundary Conditions =============
Unspecified.
=================== Results ===================
W def = 371.24
Svm max = 166.81
Evm max = 0.09 %
Ux max = 7.49e-02
Ux min = -7.49e-02
Uy max = 0.00e+00
Uy min = -9.28e-01
=================== TicTac ===================
Mesh : 235.252 ms
Boundary Conditions : 748.873 µs
Matrix : 1.942 s
Solver : 1.997 s
Display : 336.066 ms
PostProcessing : 132.556 ms
Resolution hyperelastic : 3.863 s
PyVista_Interface : 11.820 s
=================== Result ===================
err W : 0.24 %
err uy : 0.68 %
12 from EasyFEA import Display, Models, plt, np, ElemType, Simulations
13 from EasyFEA.Geoms import Domain
14
15 if __name__ == "__main__":
16 Display.Clear()
17
18 # ----------------------------------------------
19 # Configuration
20 # ----------------------------------------------
21
22 # geom
23 dim = 2
24 L = 120 # mm
25 h = 13
26 I = h**4 / 12 # mm4
27
28 # model
29 E = 210000 # MPa (Young's modulus)
30 v = 0.3 # Poisson's ratio
31 coef = 1
32
33 # load
34 load = 800 # N
35
36 # expected results
37 W_an = 2 * load**2 * L / E / h**2 * (L**2 / h**2 + (1 + v) * 3 / 5) # mJ
38 uy_an = load * L**3 / (3 * E * I)
39
40 # ----------------------------------------------
41 # Mesh
42 # ----------------------------------------------
43
44 N = 3
45 meshSize = h / N
46
47 domain = Domain((0, 0), (L, h), meshSize)
48
49 if dim == 2:
50 mesh = domain.Mesh_2D([], ElemType.QUAD9, isOrganised=True)
51 else:
52 mesh = domain.Mesh_Extrude(
53 [], [0, 0, -h], [N], ElemType.HEXA27, isOrganised=True
54 )
55
56 nodes_x0 = mesh.Nodes_Conditions(lambda x, y, z: x == 0)
57 nodes_xL = mesh.Nodes_Conditions(lambda x, y, z: x == L)
58
59 # ----------------------------------------------
60 # Simulation
61 # ----------------------------------------------
62
63 material = Models.ElasIsot(dim, E, v, planeStress=True, thickness=h)
64 simu = Simulations.ElasticSimu(mesh, material)
65
66 simu.add_dirichlet(nodes_x0, [0] * dim, simu.Get_unknowns())
67 simu.add_surfLoad(nodes_xL, [-load / h**2], ["y"])
68
69 sol = simu.Solve()
70 simu.Save_Iter()
71
72 uy_num = -simu.Result("uy").min()
73 W_num = simu._Calc_Psi_Elas()
74
75 # ----------------------------------------------
76 # Results
77 # ----------------------------------------------
78 print(simu)
79
80 Display.Section("Result")
81
82 print(f"err W : {np.abs(W_an - W_num) / W_an * 100:.2f} %")
83
84 print(f"err uy : {np.abs(uy_an - uy_num) / uy_an * 100:.2f} %")
85
86 Display.Plot_Mesh(simu, h / 2 / np.abs(sol).max())
87 Display.Plot_BoundaryConditions(simu)
88 Display.Plot_Result(simu, "uy", nodeValues=True, coef=1 / coef, ncolors=20)
89 Display.Plot_Result(simu, "Svm", plotMesh=True, ncolors=11)
90
91 plt.show()
Total running time of the script: (0 minutes 0.236 seconds)