geometricallyExactBeam2DquadraticTest.py

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

 1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 2# This is an EXUDYN example
 3#
 4# Details:  The 3-node (quadratic) ObjectBeamGeometricallyExact2D (#2208): nodes at the start, the middle
 5#           and the end of the element, in this order. (1) A cantilever with a small tip force: one
 6#           quadratic element gives the Timoshenko deflection F*L^3/(3*EI) + F*L/GA = -3.43333e-4 to 1e-7,
 7#           the size of the geometrically nonlinear part (the linear element: -3.4203e-4 with 8). (2) The same cantilever bent into a half
 8#           circle by a tip moment pi*EI/L, tip at y = 2L/pi = 0.636620: 0.616191, 0.635647, 0.636563 and
 9#           0.636616 with 1, 2, 4 and 8 quadratic elements - fourth order - against 0.707107, 0.653281 and
10#           0.640729 with 2, 4 and 8 linear elements. (3) A pendulum of 4 quadratic elements falling under
11#           gravity through MarkerBodyMass, hinged at one end: tip y -0.001862 at t = 1, 8 and 16 elements
12#           give -0.001827 and -0.001830.
13#
14# Author:   Johannes Gerstmayr
15# Date:     2026-09-30
16#
17# 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.
18#
19#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
20
21import exudyn as exu
22from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
23
24import numpy as np
25
26testIsActive = exu.sys.get('testIsActive', False)
27
28L = 1; EI = 100; GA = 1e4; EA = 1e5
29
30
31def Beam(nElements, force=0, moment=0, dynamic=False):
32    """the tip position of a beam of quadratic elements: a clamped cantilever, or a pendulum under gravity"""
33    SC = exu.SystemContainer(); mbs = SC.AddSystem()
34    nodes = [mbs.AddNode(NodeRigidBody2D(referenceCoordinates=[L*i/(2*nElements), 0, 0])) for i in range(2*nElements+1)]
35    elements = [mbs.AddObject(ObjectBeamGeometricallyExact2D(nodeNumbers=nodes[2*k:2*k+3], length=L/nElements,
36                              massPerLength=1, crossSectionInertia=0.01, bendingStiffness=EI,
37                              axialStiffness=EA, shearStiffness=GA)) for k in range(nElements)]
38    mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=mbs.AddNode(NodePointGround()), coordinate=0))
39    for i in ([0, 1] if dynamic else [0, 1, 2]):
40        mbs.AddObject(ObjectConnectorCoordinate(markerNumbers=[mGround, mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nodes[0], coordinate=i))]))
41    simulationSettings = exu.SimulationSettings()
42    simulationSettings.timeIntegration.newton.useModifiedNewton = False #Just for the test; modified Newton is usually faster
43    simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
44    if dynamic:
45        for e in elements:
46            mbs.AddLoad(LoadMassProportional(markerNumber=mbs.AddMarker(MarkerBodyMass(bodyNumber=e)), loadVector=[0, -9.81, 0]))
47        mbs.Assemble()
48        simulationSettings.timeIntegration.numberOfSteps = 1000
49        simulationSettings.timeIntegration.endTime = 1
50        simulationSettings.solution.file.write = False
51        mbs.SolveDynamic(simulationSettings)
52    else:
53        mTip = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nodes[-1]))
54        mbs.AddLoad(LoadForceVector(markerNumber=mTip, loadVector=[0, force, 0]))
55        mbs.AddLoad(LoadTorqueVector(markerNumber=mTip, loadVector=[0, 0, moment]))
56        mbs.Assemble()
57        simulationSettings.staticSolver.numberOfLoadSteps = 10 if moment != 0 else 1
58        mbs.SolveStatic(simulationSettings)
59    return mbs.GetNodeOutput(nodes[-1], exu.OutputVariableType.Position)
60
61
62force = -0.1
63tipForce = Beam(1, force=force)[1]
64exu.Print('small tip force, one element: deflection', tipForce, ', Timoshenko', force*L**3/(3*EI)+force*L/GA)
65testResult = 1000*tipForce
66
67for nElements in [2, 4]:
68    tipMoment = Beam(nElements, moment=np.pi*EI/L)[1]
69    exu.Print('half circle,', nElements, 'elements: tip y', tipMoment, ', exact', 2*L/np.pi)
70    testResult += tipMoment
71
72tipPendulum = Beam(4, dynamic=True)[1]
73exu.Print('pendulum, 4 elements: tip y at t = 1:', tipPendulum)
74testResult += 100*tipPendulum
75
76exu.Print('solution of geometricallyExactBeam2DquadraticTest=', testResult)
77exu.sys['testResult'] = testResult