rightAngleFrame.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Lateral buckling of the right-angle frame (Argyris et al. 1979; Simo and Vu-Quoc 1986;
  5#           Nachbagauer and Gerstmayr 2013), with ObjectANCFBeam or ObjectBeamGeometricallyExact:
  6#           two arms of length 0.24 m with a thin rectangular cross section 30 x 0.6 mm, clamped at the
  7#           end of the first arm, joined rigidly at the corner. The tip of the second arm is DRIVEN by
  8#           displacement: a coordinate constraint prescribes its displacement along the first arm, and
  9#           its reaction force is the load P. So the static solver passes the buckling point, where a
 10#           load-driven computation stops at the critical load 1.088 N, and follows the post-buckling
 11#           path, on which P stays close to 1.088 N while the tip moves out of the plane. A small
 12#           out-of-plane force is the imperfection.
 13#           With ObjectANCFBeam (useGeometricallyExact = False) the static solver stops at about 3 % of
 14#           the drive (P = 0.31 N): there the Newton iteration converges only from very close to the
 15#           solution, although its Jacobian is consistent and the constrained tangent stiffness is
 16#           positive definite; neither more elements, another penalty factor, another linear solver
 17#           nor more load steps change it, and without the imperfection the drive is completed (#2763).
 18#           K. Nachbagauer, J. Gerstmayr. Structural and Continuum Mechanics Approaches for a 3D Shear
 19#           Deformable ANCF Beam Finite Element: Application to Buckling and Nonlinear Dynamic Examples,
 20#           Journal of Computational and Nonlinear Dynamics, Vol. 9(1), pp. 011013-1 - 011013-8, 2013.
 21#
 22# Author:   Johannes Gerstmayr
 23# Date:     2023-05-02, displacement-driven 2026-10-01
 24#
 25# 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.
 26#
 27#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 28
 29import exudyn as exu
 30from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 31import exudyn.graphics as graphics
 32
 33import numpy as np
 34
 35testIsActive = exu.sys.get('testIsActive', False)
 36
 37useGeometricallyExact = True    #False: ObjectANCFBeam, whose static solve stops at about 3 % of the drive (#2763, see Details)
 38nElements = 16                  #per arm
 39uxMax = 2e-3                    #prescribed displacement of the tip along the first arm (x); 200 steps of 1e-5 follow the post-buckling path, much larger steps may jump to another branch
 40imperfection = 1e-4             #out-of-plane force at the tip [N]
 41numberOfLoadSteps = 200
 42
 43L = 0.24                #length of each arm
 44w = 0.0006              #thickness, out of the plane of the frame
 45h = 0.03                #height, in the plane of the frame
 46Em = 7.124e10; nu = 0.31; Gm = Em/(2*(1+nu))
 47A = h*w; Iyy = h*w**3/12; Izz = w*h**3/12
 48J = h*w**3/3            #torsion of a thin rectangle
 49rho = 1e2               #not used by the static solver
 50
 51SC = exu.SystemContainer()
 52mbs = SC.AddSystem()
 53
 54sectionData = exu.BeamSection()
 55sectionData.stiffnessMatrix = np.diag([Em*A, Gm*A, Gm*A, Gm*J, Em*Iyy, Em*Izz])
 56sectionData.inertia = rho*J*np.eye(3)
 57sectionData.massPerLength = rho*A
 58
 59sectionGeometry = exu.BeamSectionGeometry()
 60sectionGeometry.polygonalPoints = exu.Vector2DList()
 61for p in [[h, -w], [h, w], [-h, w], [-h, -w]]:
 62    sectionGeometry.polygonalPoints.Append(p)
 63
 64lElem = L/nElements
 65rotZ = RotationMatrixZ(-np.pi/2)
 66
 67def Arm(position0, direction, rotation, slopes):
 68    """the nodes and elements of one arm; returns the first and the last node"""
 69    nodes = []
 70    for k in range(nElements+1):
 71        p = list(np.array(position0) + k*lElem*np.array(direction))
 72        if useGeometricallyExact:
 73            nodes += [mbs.AddNode(NodeRigidBodyEP(referenceCoordinates=p+list(RotationMatrix2EulerParameters(rotation))))]
 74        else:
 75            nodes += [mbs.AddNode(NodePointSlope23(referenceCoordinates=p+slopes))]
 76    for k in range(nElements):
 77        if useGeometricallyExact:
 78            mbs.AddObject(ObjectBeamGeometricallyExact(nodeNumbers=[nodes[k], nodes[k+1]], length=lElem,
 79                                                       sectionData=sectionData,
 80                                                       visualization=VBeam3D(sectionGeometry=sectionGeometry)))
 81        else:
 82            mbs.AddObject(ObjectANCFBeam(nodeNumbers=[nodes[k], nodes[k+1]], length=lElem,
 83                                         sectionData=sectionData, crossSectionPenaltyFactor=[1, 1, 1],
 84                                         visualization=VANCFBeam(sectionGeometry=sectionGeometry)))
 85    return nodes[0], nodes[-1]
 86
 87nInit, n1 = Arm([0, 0, 0], [1, 0, 0], np.eye(3), [0, 1, 0, 0, 0, 1])      #first arm, along x
 88nInitB, n1B = Arm([L, 0, 0], [0, -1, 0], rotZ, [1, 0, 0, 0, 0, 1])       #second arm, along -y
 89
 90nGround = mbs.AddNode(NodePointGround(visualization=VNodePointGround(show=False)))
 91mCoordinateGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))
 92
 93#clamping at the start of the first arm, and the rigid corner
 94if useGeometricallyExact:
 95    mbs.AddObject(GenericJoint(markerNumbers=[mbs.AddMarker(MarkerNodeRigid(nodeNumber=nGround)),
 96                                              mbs.AddMarker(MarkerNodeRigid(nodeNumber=nInit))],
 97                               visualization=VGenericJoint(axesRadius=2*w, axesLength=2*w)))
 98else:
 99    for i in range(9):
100        mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround,
101                                                          mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nInit, coordinate=i))]))
102mbs.AddObject(GenericJoint(markerNumbers=[mbs.AddMarker(MarkerNodeRigid(nodeNumber=n1)),
103                                          mbs.AddMarker(MarkerNodeRigid(nodeNumber=nInitB, localHT=exu.HT(rotation=rotZ)))],
104                           visualization=VGenericJoint(axesRadius=2*w, axesLength=2*w)))
105
106#the tip is driven along x; the reaction force of the constraint is the load P
107def UFoffset(mbs, t, itemNumber, lOffset):
108    return uxMax*t
109
110oDrive = mbs.AddObject(CoordinateConstraint(markerNumbers=[mCoordinateGround,
111                                                           mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=n1B, coordinate=0))],
112                                            offsetUserFunction=UFoffset))
113mbs.AddLoad(Force(markerNumber=mbs.AddMarker(MarkerNodePosition(nodeNumber=n1B)), loadVector=[0, 0, imperfection]))
114
115sForce = mbs.AddSensor(SensorObject(objectNumber=oDrive, outputVariableType=exu.OutputVariableType.Force,
116                                    storeInternal=True, writeToFile=False))
117sDisp = mbs.AddSensor(SensorNode(nodeNumber=n1B, outputVariableType=exu.OutputVariableType.Displacement,
118                                 storeInternal=True, writeToFile=False))
119mbs.Assemble()
120
121simulationSettings = exu.SimulationSettings()
122simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
123simulationSettings.staticSolver.useLoadFactor = False   #the drive follows the user function of t in [0, 1]
124simulationSettings.staticSolver.numberOfLoadSteps = numberOfLoadSteps
125simulationSettings.staticSolver.newton.relativeTolerance = 1e-8
126simulationSettings.staticSolver.newton.absoluteTolerance = 1e-7
127simulationSettings.solution.file.write = False
128simulationSettings.solution.sensors.writePeriod = 1/numberOfLoadSteps
129
130SC.visualizationSettings.bodies.beams.axialTiling = 50
131SC.visualizationSettings.view0.scene.drawWorldBasis = True
132SC.visualizationSettings.view0.scene.worldBasisSize = 0.1
133
134if not testIsActive:
135    SC.renderer.Start()
136
137mbs.SolveStatic(simulationSettings)
138
139if not testIsActive:
140    SC.renderer.DoIdleTasks()
141    SC.renderer.Stop()
142
143force = mbs.GetSensorStoredData(sForce)[:, 1]       #the load P on the tip: the force of the drive
144uz = mbs.GetSensorStoredData(sDisp)[:, 3]           #out-of-plane displacement of the tip
145testResult = 0
146for displacement in [1e-3, 5e-3, 1e-2]:
147    forceAt = np.interp(displacement, np.abs(uz), force)
148    exu.Print('out-of-plane tip displacement', displacement, 'at P =', round(forceAt, 4))
149    testResult += forceAt
150exu.Print('at the end: P =', round(force[-1], 4), ', out-of-plane tip displacement', round(uz[-1], 5))
151testResult += force[-1] + abs(uz[-1])
152
153exu.Print('solution of rightAngleFrame=', testResult)
154exu.sys['testResult'] = testResult
155
156if not testIsActive:
157    mbs.PlotSensor(sensorNumbers=[np.column_stack((np.abs(uz), force)), np.array([[0, 1.088], [np.abs(uz).max(), 1.088]])],
158                   components=[0, 0], componentsX=[-1, -1], labels=['computed', 'critical load 1.088 N'],
159                   xLabel='out-of-plane tip displacement', yLabel='P')