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 %


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

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

Gallery generated by Sphinx-Gallery