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