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

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

Gallery generated by Sphinx-Gallery