Note
Go to the end to download the full example code.
Elas8#
A cantilever cylinder undergoing torsion, loaded either by a surface traction.
phi analytical = 1.556930e-07 rad
phi numerical = 1.556930e-07 rad
relative error = 0.000 %
13 import numpy as np
14
15 from EasyFEA import Terminal, ElemType, Models, Simulations, PyVista
16 from EasyFEA.Geoms import Circle
17
18 if __name__ == "__main__":
19 Terminal.Clear()
20
21 # --------------------------------------------
22 # Mesh
23 # --------------------------------------------
24
25 L = 100 # mm
26 D = 30
27
28 circle = Circle((0, 0), D, n=(1, 0, 0))
29 circle.meshSize = circle.length / 20
30
31 mesh = circle.Mesh_Extrude([], [L, 0, 0], [L // circle.meshSize], ElemType.PRISM6)
32
33 PyVista.Plot_Mesh(mesh).show()
34
35 # --------------------------------------------
36 # Simu
37 # --------------------------------------------
38
39 model = Models.Elastic.Isotropic(3, E=210000, v=0.3)
40 simu = Simulations.Elastic(mesh, model)
41
42 nodesX0 = mesh.Nodes_Conditions(lambda x, y, z: x == 0)
43 simu.add_dirichlet(nodesX0, [0] * 3, ["x", "y", "z"])
44
45 nodesXL = mesh.Nodes_Conditions(lambda x, y, z: x == L)
46 yC, zC = 0, 0
47 T = 10 # N.mm
48 J = np.pi * D**4 / 32 # polar second moment of area
49 p = T / J # so that the traction t = p * r integrates to a moment M_x = c
50 simu.add_surfLoad(
51 nodesXL,
52 [
53 lambda x, y, z: -p * (z - zC),
54 lambda x, y, z: p * (y - yC),
55 ],
56 ["y", "z"],
57 )
58
59 u = simu.Solve()
60
61 # --------------------------------------------
62 # Results
63 # --------------------------------------------
64 # Saint-Venant torsion of a circular shaft: phi = T * L / (G * J)
65 phi_ana = T * L / (model.get_mu() * J) # rad
66
67 # numerical twist from the tip-face nodes: a rigid rotation theta about x
68 # gives uy = -theta*dz, uz = theta*dy, so theta = (dy*uz - dz*uy) / r^2
69 uy = simu.Result("uy")[nodesXL]
70 uz = simu.Result("uz")[nodesXL]
71 dy = mesh.coord[nodesXL, 1] - yC
72 dz = mesh.coord[nodesXL, 2] - zC
73 r2 = dy**2 + dz**2
74 phi_num = np.mean((dy * uz - dz * uy) / r2, where=r2 > 2)
75
76 err = abs(phi_num - phi_ana) / abs(phi_ana) * 100
77 Terminal.MyPrint(f"\nphi analytical = {phi_ana:.6e} rad")
78 Terminal.MyPrint(f"\nphi numerical = {phi_num:.6e} rad")
79 Terminal.MyPrint(f"\nrelative error = {err:.3f} %")
80
81 PyVista.Plot_BoundaryConditions(simu).show()
82
83 pltr = PyVista._Plotter(shape=(3, 1))
84
85 PyVista.Plot(simu, "uy", plotMesh=True, plotter=pltr, verticalColobar=False)
86 pltr.subplot(1, 0)
87 PyVista.Plot(simu, "uz", plotMesh=True, plotter=pltr, verticalColobar=False)
88 pltr.subplot(2, 0)
89 PyVista.Plot(
90 simu, "displacement_norm", plotMesh=True, plotter=pltr, verticalColobar=False
91 )
92 pltr.show()
Total running time of the script: (0 minutes 0.867 seconds)