Note

Go to the end to download the full example code.

Elas6#

A bi-fixed beam undergoing bending deformation.

Elas6
  • Elas6
  • $uy^{e}$
err beam model : 0.00 %
err 3d model : 0.45 %


 12 from EasyFEA import (
 13     Display,
 14     Models,
 15     plt,
 16     np,
 17     Mesher,
 18     ElemType,
 19     Simulations,
 20     PyVista,
 21 )
 22 from EasyFEA.Geoms import Point, Points, Line
 23
 24 if __name__ == "__main__":
 25     Display.Clear()
 26
 27     # ----------------------------------------------
 28     # Configuration
 29     # ----------------------------------------------
 30
 31     # geom
 32     L = 500  # mm
 33     hx = 13
 34     hy = 20
 35     e = 2
 36
 37     # model
 38     E = 210000  # MPa (Young's modulus)
 39     v = 0.3  # Poisson's ratio
 40     coef = 1
 41
 42     # load
 43     load = 800  # N
 44
 45     # ----------------------------------------------
 46     # Section
 47     # ----------------------------------------------
 48
 49     elemType = ElemType.TETRA4
 50
 51     def DoSym(p: Point, n: np.ndarray) -> Point:
 52         pc = p.copy()
 53         pc.Symmetry(n=n)
 54         return pc
 55
 56     p1 = Point(-hx / 2, -hy / 2)
 57     p2 = Point(hx / 2, -hy / 2)
 58     p3 = Point(hx / 2, -hy / 2 + e)
 59     p4 = Point(e / 2, -hy / 2 + e, r=e)
 60     p5 = DoSym(p4, (0, 1))
 61     p6 = DoSym(p3, (0, 1))
 62     p7 = DoSym(p2, (0, 1))
 63     p8 = DoSym(p1, (0, 1))
 64     p9 = DoSym(p6, (1, 0))
 65     p10 = DoSym(p5, (1, 0))
 66     p11 = DoSym(p4, (1, 0))
 67     p12 = DoSym(p3, (1, 0))
 68     section = Points([p1, p2, p3, p4, p5, p6, p7, p8, p9, p10, p11, p12], e / 1)
 69     meshSection = Mesher().Mesh_2D(section)
 70
 71     section.Get_Contour().Plot()
 72
 73     section.Rotate(-90, direction=(0, 1))
 74
 75     # ----------------------------------------------
 76     # Mesh
 77     # ----------------------------------------------
 78
 79     mesher = Mesher()
 80     fact = mesher._factory
 81     surfaces = mesher._Surfaces(section, [], elemType)[0]
 82
 83     layers = [50]
 84     # layers = [2*L/meshSize/10]
 85
 86     mesher._Extrude(surfaces, [-L / 2, 0, 0], elemType, layers)
 87     mesher._Extrude(surfaces, [L / 2, 0, 0], elemType, layers)
 88
 89     mesher._Set_PhysicalGroups()
 90     mesher._Mesh_Generate(3, elemType)
 91     mesh = mesher._Mesh_Get_Mesh()
 92
 93     nodes_fixed = mesh.Nodes_Conditions(lambda x, y, z: (x == -L / 2) | (x == L / 2))
 94     nodes_load = mesh.Nodes_Line(Line(p7, p8))
 95
 96     # ----------------------------------------------
 97     # Beam simulation
 98     # ----------------------------------------------
 99
100     beam = Models.BeamElasIsot(
101         2, Line(Point(-L / 2), Point(L / 2), L / 10), meshSection, E, v
102     )
103
104     mesh_beam = Mesher().Mesh_Beams([beam], ElemType.SEG3)
105
106     beams = Models.BeamStructure([beam])
107     simu_beam = Simulations.BeamSimu(mesh_beam, beams)
108
109     simu_beam.add_dirichlet(
110         mesh_beam.Nodes_Conditions(lambda x, y, z: (x == -L / 2) | (x == L / 2)),
111         [0] * simu_beam.Get_dof_n(),
112         simu_beam.Get_unknowns(),
113     )
114     simu_beam.add_neumann(
115         mesh_beam.Nodes_Conditions(lambda x, y, z: (x == 0)), [-load], ["y"]
116     )
117     simu_beam.Solve()
118
119     Display.Plot_Result(simu_beam, "uy")
120
121     u_an = load * L**3 / (192 * E * beam.Iz)
122
123     uy_1d = np.abs(simu_beam.Result("uy").min())
124
125     print(f"err beam model : {np.abs(u_an - uy_1d) / u_an * 100:.2f} %")
126
127     # ----------------------------------------------
128     # Elastic simulation
129     # ----------------------------------------------
130
131     material = Models.ElasIsot(3, E, v)
132     simu = Simulations.ElasticSimu(mesh, material)
133
134     simu.add_dirichlet(nodes_fixed, [0] * 3, simu.Get_unknowns())
135     simu.add_lineLoad(nodes_load, [-load / hx], ["y"])
136     sol = simu.Solve()
137
138     uy_3d = np.abs(simu.Result("uy").min())
139
140     print(f"err 3d model : {np.abs(u_an - uy_3d) / u_an * 100:.2f} %")
141
142     # ----------------------------------------------
143     # Results
144     # ----------------------------------------------
145
146     plotter = PyVista._Plotter(shape=(2, 1))
147     PyVista.Plot(
148         simu,
149         "ux",
150         coef=1 / coef,
151         nColors=20,
152         plotMesh=True,
153         plotter=plotter,
154         verticalColobar=False,
155     )
156     PyVista.Plot_BoundaryConditions(simu, plotter=plotter)
157     plotter.subplot(1, 0)
158     PyVista.Plot(
159         simu, "uy", coef=1 / coef, nColors=20, plotter=plotter, verticalColobar=False
160     )
161     plotter.show()
162
163     plt.show()

Total running time of the script: (0 minutes 1.863 seconds)

Gallery generated by Sphinx-Gallery