Note
Go to the end to download the full example code.
LinearElasticity1#
A cantilever beam undergoing bending deformation.
Note that this simulation is also performed in examples/Elastic/Elas1.py.
14 from EasyFEA import Display, ElemType, Models, Simulations, np
15 from EasyFEA.fem import Field, BiLinearForm, FeArray, Sym_Grad, Trace
16 from EasyFEA.Geoms import Domain
17
18 if __name__ == "__main__":
19 Display.Clear()
20
21 # ----------------------------------------------
22 # Configuration
23 # ----------------------------------------------
24
25 dim = 2
26
27 # geom
28 L = 120 # mm
29 h = 13
30
31 # load
32 F = -800 # N
33
34 # model
35 E = 210000 # MPa
36 v = 0.3
37
38 elastic = Models.ElasIsot(dim, E, v, planeStress=True, thickness=h)
39 lmbda = elastic.get_lambda()
40 mu = elastic.get_mu()
41
42 # ----------------------------------------------
43 # Mesh
44 # ----------------------------------------------
45
46 contour = Domain((0, 0), (L, h), h / 3)
47
48 if dim == 2:
49 mesh = contour.Mesh_2D([], ElemType.QUAD9, isOrganised=True)
50 else:
51 mesh = contour.Mesh_Extrude(
52 [], [0, 0, h], [3], ElemType.HEXA27, isOrganised=True
53 )
54
55 nodesX0 = mesh.Nodes_Conditions(lambda x, y, z: x == 0)
56 nodesXL = mesh.Nodes_Conditions(lambda x, y, z: x == L)
57
58 # ----------------------------------------------
59 # Formulations
60 # ----------------------------------------------
61
62 field = Field(mesh.groupElem, dim)
63
64 def S(u: Field) -> FeArray:
65 Eps = Sym_Grad(u)
66 return 2 * mu * Eps + lmbda * Trace(Eps) * np.eye(dim)
67
68 @BiLinearForm
69 def ComputeK(u: Field, v: Field):
70 Sig = S(u)
71 Eps = Sym_Grad(v)
72 return Sig.ddot(Eps)
73
74 weakForms = Models.WeakForms(field, ComputeK)
75
76 # ----------------------------------------------
77 # Simulations
78 # ----------------------------------------------
79
80 simu = Simulations.WeakFormSimu(mesh, weakForms)
81
82 simu.add_dirichlet(nodesX0, [0] * dim, simu.Get_unknowns())
83 simu.add_surfLoad(nodesXL, [F / h**2], ["y"])
84
85 simu.Solve()
86
87 # ----------------------------------------------
88 # Results
89 # ----------------------------------------------
90
91 Display.Plot_Result(simu, "uy", 10, plotMesh=True)
92
93 def von_mises_stress(u):
94 Stress = S(u)
95 if dim == 2:
96 xx = Stress[..., 0, 0]
97 yy = Stress[..., 1, 1]
98 xy = Stress[..., 0, 1]
99 return np.sqrt(xx**2 + yy**2 - xx * yy + 3 * xy**2)
100 else:
101 Stress_d = Stress - 1 / 3 * Trace(Stress) * np.eye(3)
102 return np.sqrt(3 / 2 * Stress_d.ddot(Stress_d))
103
104 Svm_e = field.Evaluate_e(von_mises_stress, simu.u)
105
106 Display.Plot_Result(simu, Svm_e, plotMesh=True, ncolors=11)
107
108 Display.plt.show()
Total running time of the script: (0 minutes 0.209 seconds)