geometricallyExactBeamCurvedTest.py

You can view and download this file on Github: geometricallyExactBeamCurvedTest.py

 1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 2# This is an EXUDYN example
 3#
 4# Details:  ObjectBeamGeometricallyExact with a curved reference configuration (#1494): the 45-degree
 5#           bend of Bathe and Bolourchi (1979), a cantilever on an eighth of a circle of radius 100 in the
 6#           x-y plane, loaded at the tip by a force P out of the plane. The element is stress-free in the
 7#           reference configuration of its nodes, so the unloaded bend keeps its shape. Tip position
 8#           with 8 elements at P = 300, 450, 600: (58.60, 22.14, 40.36), (52.05, 18.40, 48.57),
 9#           (46.98, 15.59, 53.46); Simo and Vu-Quoc (1986), 8 elements: (58.84, 22.33, 40.08),
10#           (52.32, 18.62, 48.39), (47.23, 15.79, 53.37). With 32 elements at P = 600: (46.90, 15.56, 53.60).
11#
12# Author:   Johannes Gerstmayr
13# Date:     2026-09-29
14#
15# Copyright:This file is part of Exudyn. Exudyn is free software. You can redistribute it and/or modify it under the terms of the Exudyn license. See 'LICENSE.txt' for more details.
16#
17#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
18
19import exudyn as exu
20from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
21
22import numpy as np
23
24testIsActive = exu.sys.get('testIsActive', False)
25
26R = 100                 #radius of the bend
27E = 1e7; G = 0.5*E      #square cross section 1 x 1
28A = 1; I = 1/12; J = 0.141
29nElements = 8
30
31
32def Bend(force):
33    """the tip position of the bend under a tip force in z"""
34    SC = exu.SystemContainer(); mbs = SC.AddSystem()
35    nodes = []
36    for i in range(nElements+1):
37        phi = np.pi/4*i/nElements
38        nodes.append(mbs.AddNode(NodeRigidBodyEP(referenceCoordinates=[R*np.sin(phi), R*(1-np.cos(phi)), 0]
39                                                 + list(RotationMatrix2EulerParameters(RotationMatrixZ(phi))))))
40    section = exu.BeamSection()
41    section.stiffnessMatrix = np.diag([E*A, G*A, G*A, G*J, E*I, E*I])
42    section.inertia = np.diag([2*I, I, I])
43    section.massPerLength = A
44    for i in range(nElements):
45        mbs.AddObject(ObjectBeamGeometricallyExact(nodeNumbers=[nodes[i], nodes[i+1]], length=R*np.pi/4/nElements,
46                                                   sectionData=section))
47    #clamped: position and e1..e3; e0 follows from the norm constraint
48    mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=mbs.AddNode(NodePointGround()), coordinate=0))
49    for i in [0, 1, 2, 4, 5, 6]:
50        mbs.AddObject(ObjectConnectorCoordinate(markerNumbers=[mGround, mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nodes[0], coordinate=i))]))
51    mbs.AddLoad(LoadForceVector(markerNumber=mbs.AddMarker(MarkerNodePosition(nodeNumber=nodes[-1])), loadVector=[0, 0, force]))
52    mbs.Assemble()
53    simulationSettings = exu.SimulationSettings()
54    simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
55    simulationSettings.staticSolver.numberOfLoadSteps = 20
56    simulationSettings.staticSolver.newton.relativeTolerance = 1e-10
57    simulationSettings.staticSolver.newton.absoluteTolerance = 1e-6 #P = 0 has no load to measure the residual against: its round-off differs between platforms
58    mbs.SolveStatic(simulationSettings)
59    return mbs.GetNodeOutput(nodes[-1], exu.OutputVariableType.Position)
60
61
62testResult = 0
63for force in [0, 300, 450, 600]:
64    pTip = Bend(force)
65    exu.Print('P =', force, ': tip position', pTip.round(4))
66    testResult += 1e-2*sum(pTip)
67
68exu.Print('solution of geometricallyExactBeamCurvedTest=', testResult)
69exu.sys['testResult'] = testResult