geometricallyExactBeamRightAngleFrame.py

You can view and download this file on Github: geometricallyExactBeamRightAngleFrame.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) with
 5#           ObjectBeamGeometricallyExact (#1499): two arms of length 0.24 m with a thin rectangular
 6#           cross section 30 x 0.6 mm, clamped at the end of the first arm, the arms joined rigidly at
 7#           the corner, an in-plane force P at the end of the second arm, along the first arm, and a
 8#           small out-of-plane force 1e-4*P as imperfection. The frame buckles out of its plane at the
 9#           critical load 1.088 N of the literature. Measured with 16 elements per arm: the out-of-plane
10#           displacement of the tip reaches 1 mm at P = 1.081 N, 5 mm at 1.088 N and 10 mm at 1.092 N;
11#           8 and 32 elements agree to 0.003 N.
12#
13# Author:   Johannes Gerstmayr
14# Date:     2026-09-30
15#
16# 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.
17#
18#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
19
20import exudyn as exu
21from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
22
23import numpy as np
24
25testIsActive = exu.sys.get('testIsActive', False)
26
27L = 0.24                #length of each arm
28w = 0.0006              #thickness, out of the plane of the frame
29h = 0.03                #height, in the plane of the frame
30E = 7.124e10; nu = 0.31; G = E/(2*(1+nu))
31A = h*w; Iyy = h*w**3/12; Izz = w*h**3/12
32J = h*w**3/3            #torsion of a thin rectangle
33nElements = 16          #per arm
34forceMax = 1.2
35imperfection = 1e-4
36
37SC = exu.SystemContainer()
38mbs = SC.AddSystem()
39
40section = exu.BeamSection()
41section.stiffnessMatrix = np.diag([E*A, G*A, G*A, G*J, E*Iyy, E*Izz])
42section.inertia = 100*J*np.eye(3)   #not used by the static solver
43section.massPerLength = 100*A
44lElement = L/nElements
45
46#first arm along x, second arm along -y
47rotationArm2 = RotationMatrixZ(-np.pi/2)
48nodesArm1 = [mbs.AddNode(NodeRigidBodyEP(referenceCoordinates=[i*lElement, 0, 0]+list(eulerParameters0)))
49             for i in range(nElements+1)]
50nodesArm2 = [mbs.AddNode(NodeRigidBodyEP(referenceCoordinates=[L, -i*lElement, 0]+list(RotationMatrix2EulerParameters(rotationArm2))))
51             for i in range(nElements+1)]
52for nodes in [nodesArm1, nodesArm2]:
53    for i in range(nElements):
54        mbs.AddObject(ObjectBeamGeometricallyExact(nodeNumbers=[nodes[i], nodes[i+1]], length=lElement, sectionData=section))
55
56mGround = mbs.AddMarker(MarkerNodeRigid(nodeNumber=mbs.AddNode(NodePointGround())))
57mbs.AddObject(GenericJoint(markerNumbers=[mGround, mbs.AddMarker(MarkerNodeRigid(nodeNumber=nodesArm1[0]))]))
58mbs.AddObject(GenericJoint(markerNumbers=[mbs.AddMarker(MarkerNodeRigid(nodeNumber=nodesArm1[-1],
59                                                                        localHT=exu.HT(rotation=rotationArm2))),
60                                          mbs.AddMarker(MarkerNodeRigid(nodeNumber=nodesArm2[0]))]))
61
62def LoadTip(mbs, t, loadVector):
63    return [forceMax*t, 0, imperfection*forceMax*t]
64
65mbs.AddLoad(LoadForceVector(markerNumber=mbs.AddMarker(MarkerNodeRigid(nodeNumber=nodesArm2[-1])), loadVector=[0, 0, 0],
66                            loadVectorUserFunction=LoadTip))
67sTip = mbs.AddSensor(SensorNode(nodeNumber=nodesArm2[-1], outputVariableType=exu.OutputVariableType.Displacement,
68                                storeInternal=True, writeToFile=False))
69mbs.Assemble()
70
71simulationSettings = exu.SimulationSettings()
72simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
73simulationSettings.staticSolver.useLoadFactor = False   #the load follows the user function of t in [0, 1]
74simulationSettings.staticSolver.numberOfLoadSteps = 240
75simulationSettings.staticSolver.newton.relativeTolerance = 1e-8
76simulationSettings.staticSolver.newton.absoluteTolerance = 1e-7 #forces of 1 N; the default 1e-10 is below the round-off
77simulationSettings.solution.file.write = False
78simulationSettings.solution.sensors.writePeriod = 1e-4
79
80mbs.SolveStatic(simulationSettings)
81
82data = mbs.GetSensorStoredData(sTip)
83force = forceMax*data[:, 0]
84uz = data[:, 3]
85testResult = 0
86for displacement in [1e-3, 5e-3, 1e-2]:
87    forceAt = np.interp(displacement, uz, force)
88    exu.Print('out-of-plane tip displacement', displacement, 'at P =', round(forceAt, 4))
89    testResult += forceAt
90testResult += uz[-1]
91exu.Print('out-of-plane tip displacement at P =', forceMax, ':', round(uz[-1], 6))
92
93exu.Print('solution of geometricallyExactBeamRightAngleFrame=', testResult)
94exu.sys['testResult'] = testResult
95
96if not testIsActive:
97    mbs.PlotSensor(sensorNumbers=[np.column_stack((uz, force))], components=[0], componentsX=[-1],
98                   xLabel='out-of-plane tip displacement', yLabel='P')