geometricallyExactBeamTest.py

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

  1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  2# This is an EXUDYN example
  3#
  4# Details:  Test models for GeometricallyExactBeam (2-node shear deformable beam, 
  5#           Lie group formulation for work of elastic forces); 
  6#           test models: cantilever beam with tip force and torque
  7#
  8# Author:   Johannes Gerstmayr
  9# Date:     2023-04-05
 10#
 11# 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.
 12#
 13#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 14
 15import exudyn as exu
 16from exudyn.utilities import * #includes itemInterface and rigidBodyUtilities
 17from math import sin, cos
 18import exudyn.graphics as graphics
 19
 20import numpy as np
 21
 22testIsActive = exu.sys.get('testIsActive', False)
 23
 24SC = exu.SystemContainer()
 25mbs = SC.AddSystem()
 26
 27compute2D = False
 28compute3D = True
 29
 30#test examples
 31#2011 MUBO, Nachbagauer Pechstein Irschik Gerstmayr (2D)
 32#2013 CND, Nachbagauer Gruber Gerstmayr (static, 3D); "Structural and Continuum Mechanics Approaches for a 3D Shear Deformable ANCF Beam Finite Element: Application to Static and Linearized Dynamic Examples"
 33### not yet: 2013 CND, Nachbagauer Gerstmayr (dynamic, 3D)
 34cases = ['CantileverLinear2011', 'Cantilever2011', 'GeneralBending2013', 'PrincetonBeamF2', 'PrincetonBeamF3', 'Eigenmodes2013']
 35nElementsList = [1,2,4,8,16,32,64,128,256,512,1024]
 36# nElementsList = [8,32, 128]
 37nElements = 8
 38
 39betaList = [0,15,30,45,60,75,90]
 40betaDegree = 45
 41
 42caseList = [0,1,2,3,4] #case 0 not working for Geometrically exact beam
 43#case=2
 44
 45testIsActive = True
 46verbose = 1*0
 47
 48useEP = True #for geometrically exact beam node
 49if useEP:
 50    NodeClass = NodeRigidBodyEP
 51    initialRotationsGE = eulerParameters0
 52else: #does not work for static case, as static solver currently (2023-04) cannot solve for Lie group nodes 
 53    NodeClass = NodeRigidBodyRotVecLG
 54    initialRotationsGE = [0,0,0]
 55
 56
 57bodyFixedLoad = False
 58testErrorSum = 0
 59case=4
 60
 61printCase = True
 62#for nElements in nElementsList:
 63#for betaDegree in betaList:
 64for case in caseList:
 65#if True:
 66    if printCase:
 67        printCase=False
 68        exu.Print('case=', case, cases[case])
 69    mbs.Reset()
 70    
 71    computeEigenmodes = False
 72    csFact = 1
 73    sectionData = exu.BeamSection()
 74    fTip = 0
 75    MxTip = 0
 76    MyTip = 0
 77
 78    ks1=1 #shear correction, torsion
 79    ks2=1 #shear correction, bending
 80    ks3=1 #shear correction, bending
 81    ff=1 #drawing factor
 82
 83    if case == 0 or case == 1:
 84        caseName = cases[case]
 85
 86        L = 2 #length of beam
 87        w = 0.1 #width of beam
 88        h = 0.5 #height Y
 89
 90        fTip = 5e5*h**3
 91        if case == 1:
 92            fTip *= 1000
 93
 94        Em = 2.07e11
 95        rho = 1e2
 96    
 97        A=h*w
 98        nu = 0.3              # Poisson ratio
 99        ks2= 10*(1+nu)/(12+11*nu)
100        ks3=ks2
101
102    elif case == 2:
103        L = 2 #length of beam
104        h = 0.2 #height Y
105        w = 0.4 #width Z of beam
106        Em = 2.07e11
107        rho = 1e2
108
109        A=h*w
110    
111        nu = 0.3              # Poisson ratio
112        ks1= 0.5768 #torsion correction factor if J=Jyy+Jzz
113        ks2= 0.8331
114        ks3= 0.7961
115    
116        MxTip = 0.5e6
117        MyTip = 2e6
118
119        csFact = 10
120    elif case == 3 or case == 4: #Princeton beam example
121        L = 0.508       #length of beam
122        h = 12.3777e-3  #height Y; 12.3777e-3 with Obrezkov's paper
123        w = 3.2024e-3   #width Z of beam
124        Em = 71.7e9
125        ks1=0.198
126        nu = 0.31
127
128        ks2=1
129        ks3=1
130        # ks2=0.9
131        # ks3=0.9
132
133
134        rho = 1e2       #unused
135        A=h*w
136    
137        MxTip = 0
138        MyTip = 0
139        if case == 3:
140            fTip = 8.896    #F2
141        elif case == 4:
142            fTip = 13.345 #F3
143        #if kk==0: exu.Print('load=', fTip)
144        
145        beta = betaDegree/180*pi #beta=0 => negative y-axis
146        bodyFixedLoad = False
147
148        csFact = 10
149    
150    Gm = Em/(2*(1+nu))      # Shear modulus
151
152    # Cross-section properties
153    Iyy = h*w**3/12 # Second moment of area of the beam cross-section
154    Izz = w*h**3/12 # Second moment of area of the beam cross-section
155    J = (Iyy+Izz)   # approximation; Polar moment of area of the beam cross-section
156
157    sectionData.stiffnessMatrix = np.diag([Em*A, Gm*A*ks2, Gm*A*ks3, Gm*J*ks1, Em*Iyy, Em*Izz])
158
159
160    rhoA = rho*A
161
162    if False:
163        #linear solution:
164        uzTip = fTip*L**3/(3*Em*Iyy)
165        exu.Print('uz linear=',uzTip)
166        uyTip = fTip*L**3/(3*Em*Izz)
167        exu.Print('uy linear=',uyTip)
168
169    sectionData.inertia= rho*J*np.eye(3)
170    sectionData.massPerLength = rhoA
171
172    sectionGeometry = exu.BeamSectionGeometry()
173
174    #points, in positive rotation sense viewing in x-direction, points in [Y,Z]-plane
175    #points do not need to be closed!
176    lp = exu.Vector2DList()
177    if True:
178        lp.Append([h*ff,-w*ff])
179        lp.Append([h*ff,w*ff])
180        lp.Append([-h*ff,w*ff])
181        lp.Append([-h*ff,-w*ff])
182
183    sectionGeometry.polygonalPoints = lp
184    #exu.Print('HERE\n',sectionGeometry.polygonalPoints)
185    nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
186    mnGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))
187 
188
189    eY=[0,1,0]
190    eZ=[0,0,1]
191    lElem = L/nElements
192    useGeometricallyExact = True
193    if compute3D:
194        if useGeometricallyExact:
195            n0 = mbs.AddNode(NodeClass(referenceCoordinates=[0,0,0]+initialRotationsGE))
196        else:
197            initialRotations = eY+eZ
198            n0 = mbs.AddNode(NodePointSlope23(referenceCoordinates=[0,0,0]+initialRotations))
199        nInit = n0
200        for k in range(nElements):
201            if useGeometricallyExact:
202                n1 = mbs.AddNode(NodeClass(referenceCoordinates=[lElem*(k+1),0,0]+initialRotationsGE))
203            
204                oBeam = mbs.AddObject(ObjectBeamGeometricallyExact(nodeNumbers=[n0,n1], length = lElem, 
205                                                                     sectionData = sectionData,
206                                                                     visualization=VBeam3D(sectionGeometry=sectionGeometry)))
207            else:
208                n1 = mbs.AddNode(NodePointSlope23(referenceCoordinates=[lElem*(k+1),0,0]+initialRotations))
209                oBeam = mbs.AddObject(ObjectANCFBeam(nodeNumbers=[n0,n1], length = lElem, 
210                                                       #testBeamRectangularSize = [h,w],
211                                                       sectionData = sectionData,
212                                                       crossSectionPenaltyFactor = [csFact,csFact,csFact],
213                                                       visualization=VANCFBeam(sectionGeometry=sectionGeometry)))
214            n0 = n1
215
216
217        mTip = mbs.AddMarker(MarkerNodeRigid(nodeNumber = n1))
218        if fTip != 0:
219            if case < 3:
220                mbs.AddLoad(Force(markerNumber=mTip, loadVector = [0,fTip,0], bodyFixed = bodyFixedLoad))
221            elif case >= 3:
222                mbs.AddLoad(Force(markerNumber=mTip, loadVector = [0,-fTip*cos(beta),fTip*sin(beta)], bodyFixed = bodyFixedLoad))
223
224        if MxTip != 0 or MyTip != 0:
225            mbs.AddLoad(Torque(markerNumber=mTip, loadVector = [MxTip, MyTip,0]))#, bodyFixed = True ))
226
227        if useGeometricallyExact:
228            nm0 = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nInit))
229            nmGround = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nGround))
230            mbs.AddObject(GenericJoint(markerNumbers=[nmGround, nm0]))
231        else:
232            for i in range(9):
233                #if i != 4 and i != 8: #exclude constraining the slope lengths
234                if True:
235                    nm0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nInit, coordinate=i))
236                    mbs.AddObject(CoordinateConstraint(markerNumbers=[mnGround, nm0]))
237
238
239    # exu.Print(mbs)
240    mbs.Assemble()
241
242    tEnd = 100     #end time of simulation
243    stepSize = 0.5*0.01*0.1    #step size; leads to 1000 steps
244
245    simulationSettings = exu.SimulationSettings()
246    simulationSettings.solution.file.writePeriod = 2e-2  #output interval general
247    simulationSettings.solution.sensors.writePeriod = 1e-1  #output interval of sensors
248    simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize) #must be integer
249    simulationSettings.timeIntegration.endTime = tEnd
250    #simulationSettings.solution.file.information = "This is the info\nNew line\n and another new line \n"
251    simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
252    #simulationSettings.timeIntegration.realtime.active=True
253    #simulationSettings.timeIntegration.realtime.factor=0.1
254
255    simulationSettings.timeIntegration.verboseMode = verbose
256    simulationSettings.staticSolver.verboseMode = verbose
257
258    simulationSettings.timeIntegration.newton.useModifiedNewton = True
259    #simulationSettings.timeIntegration.newton.numericalDifferentiation.minimumCoordinateSize = 1e0
260
261    #simulationSettings.timeIntegration.newton.numericalDifferentiation.relativeEpsilon = 1e-4
262    simulationSettings.timeIntegration.newton.relativeTolerance = 1e-6
263
264    # simulationSettings.show.computationTime = True
265    simulationSettings.linearSolver.solverType = exu.LinearSolverType.EigenSparse
266
267    simulationSettings.staticSolver.numberOfLoadSteps = 5
268    simulationSettings.staticSolver.adaptiveStep = True
269    #simulationSettings.staticSolver.stabilizerODE2term = 100
270
271    if useGeometricallyExact:
272        # simulationSettings.staticSolver.newton.numericalDifferentiation.forODE2 = True
273        # simulationSettings.staticSolver.newton.numericalDifferentiation.relativeEpsilon = 1e-5
274        if case == 0:
275            simulationSettings.staticSolver.newton.relativeTolerance = 1e-4
276            simulationSettings.staticSolver.newton.absoluteTolerance = 1e-5
277            simulationSettings.staticSolver.numberOfLoadSteps = 1 #otherwise makes problems
278
279    if nElements > 32 and case==0: #change tolerance, because otherwise no convergence
280        simulationSettings.staticSolver.newton.relativeTolerance = 1e-6
281    if case == 1: #tolerance changed from 1e-8 to 5e-10 to achieve values of paper (1024 has difference at last digit in paper)
282        simulationSettings.staticSolver.newton.relativeTolerance = 0.5e-9
283
284
285
286    #add some drawing parameters for this example
287    SC.visualizationSettings.nodes.drawNodesAsPoint=False
288    SC.visualizationSettings.nodes.defaultSize=0.01
289
290    SC.visualizationSettings.bodies.beams.axialTiling = 50
291    SC.visualizationSettings.view0.scene.drawWorldBasis = True
292    SC.visualizationSettings.view0.scene.worldBasisSize = 0.1
293    SC.visualizationSettings.openGL.multiSampling = 4
294
295
296    # [M, K, D] = exu.solver.ComputeLinearizedSystem(mbs, simulationSettings, useSparseSolver=True)
297    # exu.Print('M=',M.round(1))
298
299    if not testIsActive:
300        SC.renderer.Start()
301        SC.renderer.DoIdleTasks()
302
303    mbs.SolveStatic(simulationSettings)
304    # mbs.SolveDynamic(simulationSettings)
305    #mbs.SolveDynamic(simulationSettings, solverType = exu.DynamicSolverType.RK44)
306
307    #check jacobian
308    if False:
309    #%%+++++++++++++++++++++++++++++++++++    
310        solver=mbs.sys['staticSolver']
311        solver.InitializeSolver(mbs, simulationSettings)
312        solver.ComputeJacobianODE2RHS(mbs)
313        J=solver.GetSystemJacobian()
314        exu.Print((1e-6*J[:14,:7]).round(3))
315        exu.Print((1e-6*J[:14,7:14]).round(3))
316
317    #%%+++++++++++++++++++++++++++++++++++    
318    if not testIsActive:
319        SC.renderer.DoIdleTasks()
320        SC.renderer.Stop() #safely close rendering window!
321    
322    ##evaluate final (=current) output values
323    uTip = mbs.GetNodeOutput(n1, exu.OutputVariableType.Displacement)
324    
325    errorFact = 1
326    if case != 1:
327        errorFact *= 100
328
329    testErrorSum += np.linalg.norm(uTip)
330
331    
332    
333    if case < 2:
334        pTip = mbs.GetNodeOutput(n1, exu.OutputVariableType.Position)
335        exu.Print('ne=',nElements, ', ux=',L-pTip[0], ', uy=',pTip[1])
336    elif case == 2:
337        rotTip = mbs.GetNodeOutput(n1, exu.OutputVariableType.Rotation)
338        exu.Print('ne=',nElements, ', u=',list(uTip))
339        # exu.Print('ne=',nElements, ', rot=',rotTip)
340    elif case == 3 or case == 4:
341        exu.Print('ne=', nElements, ', beta=', round(beta*180/pi,1), ', u=',uTip.round(7))
342
343
344exu.Print('Solution of geometricallyExactBeamTest=', testErrorSum)
345exu.sys['testResult'] = testErrorSum
346
347
348#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
349#all results are taken from ANCFBeam (shear deformable 2-node 3D beam):
350#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
351
352# case= 0/CantileverLinear2011 
353#NachbagauerPechsteinIrschikGerstmayrMUBO2011 (2D):
354# ne=1,   9.12273046e–8, 6.16666566e–4, 0.000193
355# ne=2,   1.61293091e–7, 7.61594059e–4, 4.831e–5
356# ne=4,   1.81763233e–7, 7.97825954e–4, 1.208e–5
357# ne=256, 1.88847418e–7, 8.09900305e–4, 2.945e–9
358#Exudyn: ksFact=1
359# ne= 1 , ux= 9.122730637578513e-08 , uy= 0.0006166665660910789 
360# ne= 2 , ux= 1.612930911054633e-07 , uy= 0.0007615940599560586 
361# ne= 4 , ux= 1.8176323512975046e-07 , uy= 0.0007978259537503566 
362# ne= 8 , ux= 1.8706537496804287e-07 , uy= 0.0008068839288072378 
363# ne= 16 , ux= 1.8840244964124508e-07 , uy= 0.0008091484226773518 
364# ne= 32 , ux= 1.887374359021976e-07 , uy= 0.0008097145461515286 
365# ne= 64 , ux= 1.888212299849812e-07 , uy= 0.0008098560770202866 
366# ne= 128 , ux= 1.8884218011550047e-07 , uy= 0.000809891459736643 
367# ne= 256 , ux= 1.8884741770364144e-07 , uy= 0.0008099003054122335
368
369
370# case= 1/Cantilever2011 
371#NachbagauerPechsteinIrschikGerstmayrMUBO2011 (2D):
372# ne=1,    0.07140274, 0.54225823, 0.168310
373# ne=2,    0.12379212, 0.65687111, 0.053697
374# ne=4,    0.14346767, 0.69593561, 0.014633
375# ne=1024, 0.15097103, 0.71056837, 2.280e–7
376
377#Exudyn: ksFact=1
378# ne= 1 , ux= 0.07140273975041422 , uy= 0.5422582285449739 
379# ne= 2 , ux= 0.12379212054619537 , uy= 0.6568711099777776 
380# ne= 4 , ux= 0.14346766617229956 , uy= 0.695935613449867 
381# ne= 8 , ux= 0.14904162148449163 , uy= 0.7068152604035266 
382# ne= 16 , ux= 0.15048521526298897 , uy= 0.709623891842095 
383# ne= 32 , ux= 0.15084943688011565 , uy= 0.7103320154655514 
384# ne= 64 , ux= 0.15094070328691145 , uy= 0.7105094267817303 
385# ne= 128 , ux= 0.15096353326024237 , uy= 0.7105538037895819 
386# ne= 256 , ux= 0.15096924149743085 , uy= 0.7105648993600513 
387# ne= 512 , ux= 0.15097066651939461 , uy= 0.7105676689547459 
388# ne= 1024 , ux= 0.15097102364723924 , uy= 0.7105683631862169 
389
390# case = 2:
391#2013 CND, Nachbagauer Gruber Gerstmayr (static, 3D); "Structural and Continuum Mechanics Approaches for a 3D Shear Deformable ANCF Beam Finite Element: Application to Static and Linearized Dynamic Examples"
392#Table 4:
393# SMF
394# 8,  1.0943e-4, 1.8638e-4, 1.8117e-2
395# 32, 1.0943e-4, 1.8625e-4, 1.8117e-2
396# ANSYS
397# 40, 1.0939e-4, 1.8646e-4, 1.8117e-2    
398#Exudyn, ksFact=10:
399# ne= 8 , u= [-0.00010900977088157404, -0.0001902100873246334, -0.01811732779800177] 
400# ne= 32 , u= [-0.00010941122286522997, -0.00018667435478355072, -0.01811739809277171] 
401# ne= 128 , u= [-0.00010943631319815239, -0.000186451835025629, -0.018117402461210096] 
402#==> in 2013 paper, element performed slightly better, especially in ux and uy terms
403
404# case = 3:
405#Princeton beam with ANSYS (Leonid Obrezkov / Aki Mikkola / Marko Matikainen et al.,
406#       Performance review of locking alleviation methods for continuum ANCF beam elements,
407#       Nonlinear Dynamics, Vol. 109, pp. 31–546, May 2022
408# beta=[0 15 30 45 60 75 90];
409if (case==3 or case == 4) and False:
410    # F2=8.896
411    # % ANSYS beam (10-199 el)
412    ANSYSF2y=np.array([1.071417630E-002,  1.061328706E-002, 1.011169630E-002,  8.837226265E-003, 6.604665004E-003, 3.538889001E-003, 0])
413    ANSYSF2z=np.array([0, 4.208232124E-002, 7.939482948E-002, 0.108987937,  0.129887616, 0.142194370, 0.146245978])
414    exu.Print('refsol ANSYS F2=8.896:\n',ANSYSF2y.round(6), '\n', ANSYSF2z.round(6))
415    # % ANSYS solid (el) (4x12x500) - finer mesh doesn't have much influence see in Size effect file
416    # ANSYS_solid_y=[1.069752828E-002 1.057180106E-002 9.938278402E-003 8.686786771E-003 6.500006282E-003 3.481999513E-003 0];
417    # ANSYS_solid_z=[0 4.101165651E-002 7.696749069E-002 0.105976311 0.127251299 0.139594740 0.143848652];
418    
419    # F3=13.345 
420    # % ANSYS beam (10-199 el)
421    ANSYSF3y=np.array([1.606423724E-002, 1.645825752E-002, 1.665873206E-002, 1.518618440E-002, 1.157837500E-002, 6.248967384E-003, 0])
422    ANSYSF3z=np.array([0,                6.435812858E-002, 0.117735994,      0.156467239,      0.181861627,      0.196097131,      0.200677707])
423    # % ANSYS solid (el) (4x12x500) - finer mesh see in Size effect file
424    #ANSYS_solid_y=[1.603700622E-002 1.637026068E-002 1.640440775E-002 1.485055210E-002 1.127173264E-002 6.062461977E-003  0])
425    #ANSYS_solid_z=[0 6.270699533E-002 0.113752002 0.153554457 0.179978534  0.192972233 0.197669499])
426    exu.Print('refsol ANSYS F3=13.345:\n',ANSYSF3y.round(6), '\n', ANSYSF3z.round(6))
427#Exudyn results for Princeton beam:
428#not exactly the same, but around the previous values with HOTINT
429#using 16 elements, csFact=10 (no influence)
430# F2=8.896
431# case= 3, PrincetonBeam
432# ne= 16 , beta= 0.0 , u= [-0.0001352 -0.0107023  0.       ]
433# ne= 16 , beta= 15.0 , u= [-0.0022414 -0.0106295  0.0421374]
434# ne= 16 , beta= 30.0 , u= [-0.0076567 -0.0101861  0.0794434]
435# ne= 16 , beta= 45.0 , u= [-0.0143664 -0.0089529  0.1089703]
436# ne= 16 , beta= 60.0 , u= [-0.0204225 -0.0067182  0.1297877]
437# ne= 16 , beta= 75.0 , u= [-0.0245093 -0.0036079  0.1420319]
438# ne= 16 , beta= 90.0 , u= [-0.0259403 -0.         0.1460608]
439
440# F3=13.345
441# case= 4, PrincetonBeam
442# ne= 16 , beta= 0.0 , u= [-0.0003039 -0.0160454  0.       ]
443# ne= 16 , beta= 15.0 , u= [-0.005319  -0.0165469  0.064622 ]
444# ne= 16 , beta= 30.0 , u= [-0.0171901 -0.0169316  0.1179818]
445# ne= 16 , beta= 45.0 , u= [-0.0303357 -0.0155488  0.1565214]
446# ne= 16 , beta= 60.0 , u= [-0.0411035 -0.0118996  0.1817173]
447# ne= 16 , beta= 75.0 , u= [-0.0479101 -0.0064334  0.1958343]
448# ne= 16 , beta= 90.0 , u= [-0.0502184 -0.         0.2003738]