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

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

Gallery generated by Sphinx-Gallery