Beam5

Frame with six beams

• 
• 
• 
• 
• 
• 
• 
• 
• 
• 

The beam's vertical axis has been selected incorrectly (collinear with the beam x-axis).
Axis [-1.  0.  0.] has been assigned for beam7.

Node 0 at [200.   0.   0.]
  ux=-1.30e+00 mm, uy=-5.51e+00 mm, rz=-4.11e-02 rad
  fx=3.87e+01 N, fy=-3.88e+01 N, cz=-9.11e+00 N.mm

Node 1 at [100. 100.   0.]
  ux=4.26e-01 mm, uy=-2.54e+00 mm, rz=-2.01e-03 rad
  fx=-3.87e+01 N, fy=3.88e+01 N, cz=1.69e+01 N.mm

Node 2 at [200.   0.   0.]
  ux=-1.30e+00 mm, uy=-5.51e+00 mm, rz=-4.11e-02 rad
  fx=-3.87e+01 N, fy=-4.72e-01 N, cz=9.11e+00 N.mm

Node 3 at [100.   0.   0.]
  ux=-8.64e-01 mm, uy=-2.10e+00 mm, rz=-1.03e-02 rad
  fx=3.87e+01 N, fy=4.72e-01 N, cz=3.81e+01 N.mm

Node 4 at [100.   0.   0.]
  ux=-8.64e-01 mm, uy=-2.10e+00 mm, rz=-1.03e-02 rad
  fx=-1.08e+00 N, fy=3.84e+01 N, cz=-5.78e+01 N.mm

import matplotlib.pyplot as plt
import numpy as np

from EasyFEA import Display, Models, Mesher, ElemType, Simulations
from EasyFEA.Geoms import Domain, Line

if __name__ == "__main__":
    Display.Clear()

    # ----------------------------------------------
    # Configuration
    # ----------------------------------------------

    # geom
    L = 100  # mm

    # model
    E = 276  # MPa
    v = 0.3
    beamDim = 3  # must be >= 2

    # ----------------------------------------------
    # Mesh
    # ----------------------------------------------

    pA = (2 * L, 0)
    pB = (L, 0)
    pC = (L, L)
    pD = (0, 0)
    pE = (0, L)

    line1 = Line(pA, pC)
    line2 = Line(pA, pB)
    line3 = Line(pB, pC)
    line4 = Line(pC, pE)
    line5 = Line(pB, pD)
    line6 = Line(pB, pE)
    lines = [line1, line2, line3, line4, line5, line6]

    contour = Domain((-4 / 2, -8 / 2), (4 / 2, 8 / 2))
    section = Mesher().Mesh_2D(contour)

    beams = [Models.Beam.Isotropic(beamDim, line, section, E, v) for line in lines]
    structure = Models.Beam.BeamStructure(beams)

    mesh = Mesher().Mesh_Beams(beams, ElemType.SEG2)
    # Display.Plot_Mesh(mesh)
    # Display.Plot_Tags(mesh)

    # ----------------------------------------------
    # Simulation
    # ----------------------------------------------
    simu = Simulations.Beam(mesh, structure)

    nodesRigi = mesh.Nodes_Point(pE)
    nodesRigi = np.append(nodesRigi, mesh.Nodes_Point(pD))
    nodesA = mesh.Nodes_Point(pA)

    # link beams at specified points
    for point in [pA, pB, pC]:
        nodes = mesh.Nodes_Point(point)
        firstNodes = nodes[0]
        others = nodes[1:]
        [simu.add_connection_hinged([firstNodes, n]) for n in others]

    simu.add_dirichlet(nodesRigi, [0, 0], ["x", "y"])
    simu.add_neumann(nodesA, [-4 * 9.81], ["y"])

    simu.Solve()

    # ----------------------------------------------
    # Results
    # ----------------------------------------------
    matrixDep = simu.Results_displacement_matrix()
    depMax = np.max(np.linalg.norm(matrixDep, axis=1))

    Display.Plot_Mesh(simu, deformFactor=10 / depMax)
    Display.Plot_Mesh(section, title="Cross section")
    Display.Plot_BoundaryConditions(simu)
    Display.Plot_Result(simu, "ux", deformFactor=5 / depMax)
    Display.Plot_Result(simu, "uy", deformFactor=5 / depMax)
    Display.Plot_Result(simu, "rz", deformFactor=5 / depMax)
    Display.Plot_Result(simu, "fx", deformFactor=5 / depMax)
    Display.Plot_Result(simu, "fy", deformFactor=5 / depMax)

    Epsilon_e_pg = simu._Calc_Epsilon_e_pg(simu.displacement)
    Internal_e = simu._Calc_InternalForces_e_pg(Epsilon_e_pg).mean(1)
    Sigma_e = simu._Calc_Sigma_e_pg(Epsilon_e_pg).mean(1)
    Display.Plot_Result(simu, Sigma_e[:, 0], title="Sxx")
    Display.Plot_Result(simu, Internal_e[:, 0], title="N")

    ux, uy, rz = simu.Result("ux"), simu.Result("uy"), simu.Result("rz")
    fx, fy, cz = simu.Result("fx"), simu.Result("fy"), simu.Result("cz")

    for i in range(5):
        print(f"\nNode {i} at {simu.mesh.coord[i]}")
        print(f"  ux={ux[i]:.2e} mm, uy={uy[i]:.2e} mm, rz={rz[i]:.2e} rad")
        print(f"  fx={fx[i]:.2e} N, fy={fy[i]:.2e} N, cz={cz[i]:.2e} N.mm")

    plt.show()