NGsolveFFRFSlidingJoint.py
You can view and download this file on Github: NGsolveFFRFSlidingJoint.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Test for sliding joint and beam attached to FEM mesh
5#
6# Author: Johannes Gerstmayr
7# Date: 2026-02-02
8#
9# 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.
10#
11#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
12
13
14import exudyn as exu
15from exudyn.utilities import *
16import exudyn.graphics as graphics
17from exudyn.FEM import HCBstaticModeSelection, FEMinterface, ObjectFFRFreducedOrderInterface, KirchhoffMaterial
18from exudyn.beams import GenerateStraightLineANCFCable
19
20SC = exu.SystemContainer()
21mbs = SC.AddSystem()
22
23import numpy as np
24
25import time
26
27import netgen.occ as occ
28from ngsolve import Mesh, Draw
29
30
31
32gravity=[0,0,-9.81]
33width = 3
34height = 1.5
35length = 4
36thickness = 0.1
37arc = 0.5
38bDim = 0.08
39bDimX = 2*bDim
40rCable = bDim*0.6
41
42carrierWidth = 0.5
43carrierHeight = 0.1
44carrierLength = width-2*thickness-0.4*bDim
45carrierT = 0.05*carrierWidth
46
47
48maxh=2*thickness #*1 gives fine solution, but takes 5 minutes to simulate
49maxhC=0.5*thickness
50curvaturesafety=1.5
51endTime=3
52stepSize = 2e-3
53
54
55rho = 2800
56Emodulus = 8e11
57nu = 0.3
58materialAlu = KirchhoffMaterial(Emodulus, nu, rho)
59
60interfaceNameList = ['ground']
61
62if True:
63 #%%++++++++++++++++++++++++++++++++++++++++++++++++
64 #extrusion geometry for frame:
65 wp = occ.WorkPlane(occ.Axes(p=(0,0,0), n=occ.X, h=occ.Y))
66 #wp.MoveTo(0,0).Line(0.5*width*2).Rotate(90).Line(height).Rotate(45).Line(0.5).Close()
67 wp.MoveTo(0,0).Line(0.5*width-arc).Arc(arc,90).Line(height-arc).Rotate(90).Line(thickness).\
68 Rotate(90).Line(height-arc).Arc(arc-thickness,-90).Line(width-2*arc).Arc(arc-thickness,-90).\
69 Line(height-arc).Rotate(90).Line(thickness).Rotate(90).Line(height-arc).\
70 Arc(arc,90).Close()
71
72 frame = wp.Face().Extrude(length)
73 frame.faces.Min((0,0,1)).name='ground'
74
75 vBox = [0.5*bDimX,0.5*bDim,0.5*bDim]
76
77 boxList = []
78 for fy in [-1,1]:
79 for fx in [0,0.5,1]:
80 name = 'box_y'+str(fy)+'_x'+str(fx)
81 interfaceNameList.append(name)
82 pBox = np.array((0.5*bDimX+(length-bDimX)*fx,
83 fy*(-0.5*width+thickness+0.5*bDim),
84 0.5*height))
85 box = occ.Box(tuple(pBox-vBox),
86 tuple(pBox+vBox))
87 if fy > 0:
88 box.faces.Max((0, 1, 0)).name = name
89 else:
90 box.faces.Min((0, 1, 0)).name = name
91 boxList.append(box)
92
93 geo = occ.OCCGeometry(frame+boxList[0]+boxList[1]+boxList[2]+boxList[3]+boxList[4]+boxList[5])
94 print('meshing ...')
95
96 if False:
97 import netgen.gui #this starts netgen gui; Press button "Visual" and activate "Auto-redraw after (sec)"; Then select "Mesh"
98
99
100 geoMesh = geo.GenerateMesh(maxh=maxh,
101 curvaturesafety=curvaturesafety,
102 )
103
104 mesh = Mesh(geoMesh)
105
106 femInterface = FEMinterface()
107 [bfM, bfK, fes] = femInterface.ImportMeshFromNGsolve(mesh,
108 density=rho, youngsModulus=Emodulus, poissonsRatio=nu,
109 boundaryNamesList=interfaceNameList,
110 meshOrder=2
111 )
112 print('nNodes=', femInterface.NumberOfNodes())
113
114 [boundaryNodesList, boundaryWeightsList] = femInterface.GetBoundaryNodeSetsAsLists()
115 femInterface.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryNodesList,
116 nEigenModes=8,
117 excludeRigidBodyMotion=True,
118 boundaryNodesWeights=boundaryWeightsList,
119 computationMode=HCBstaticModeSelection.RBE2)
120
121 print('eigenfrequencies frame (Hz):\n',np.round(femInterface.GetEigenFrequenciesHz(),2),sep='')
122 femInterface.ComputePostProcessingModesNGsolve(fes, materialAlu)
123
124
125 createFFRFObjectDict0 = mbs.CreateFFRFReducedOrderObject(name='frame',
126 femInterface=femInterface,
127 stiffnessProportionalDamping=5e-4,
128 gravity=gravity)
129 mFrameGround = createFFRFObjectDict0['frame:ground']
130 mFrameBoxL = []
131 mFrameBoxR = []
132 mFrameBoxL.append(createFFRFObjectDict0['frame:box_y-1_x0'])
133 mFrameBoxL.append(createFFRFObjectDict0['frame:box_y-1_x0.5'])
134 mFrameBoxL.append(createFFRFObjectDict0['frame:box_y-1_x1'])
135 mFrameBoxR.append(createFFRFObjectDict0['frame:box_y1_x0'])
136 mFrameBoxR.append(createFFRFObjectDict0['frame:box_y1_x0.5'])
137 mFrameBoxR.append(createFFRFObjectDict0['frame:box_y1_x1'])
138
139 mFrameBox = [mFrameBoxL, mFrameBoxR]
140
141 #%%++++++++++++++++++++++++++++++++++++++++++++++++
142
143 gFloor = graphics.CheckerBoard(point=[0,0,0.],size=8)
144 oGround = mbs.CreateGround(graphicsDataList=[gFloor])
145
146 mbs.CreateGenericJoint(bodyNumbers=[mFrameGround,oGround])
147
148#%%++++++++++++++++++++++++++++++++++++++++++++++++
149interfaceNameListC = ['attachment']
150
151wp = occ.WorkPlane(occ.Axes(p=(0,0,0), n=-occ.Y, h=occ.X))
152wp.MoveTo(0,0).LineTo(carrierT,0)\
153 .LineTo(carrierT,carrierHeight-carrierT)\
154 .LineTo(carrierWidth-carrierT,carrierHeight-carrierT)\
155 .LineTo(carrierWidth-carrierT,0)\
156 .LineTo(carrierWidth,0)\
157 .LineTo(carrierWidth,carrierHeight)\
158 .LineTo(0 ,carrierHeight)\
159 .Close()
160
161carrier = wp.Face().Extrude(carrierLength).Move((-0.5*carrierWidth,0.5*carrierLength,0))
162#frame.faces.Min((0,0,1)).name='ground'
163# carrier.faces.Max((0,0,1)).name='attachment'
164# carrier.faces.Min((0,0,1)).name='base'
165
166pFront = np.array([0.,-carrierLength*0.5+carrierT*0.5,carrierHeight*0.5])
167fBox = np.array([0.5*carrierWidth,0.5*carrierT,0.5*carrierHeight])
168frontL = occ.Box(tuple(pFront-fBox),
169 tuple(pFront+fBox))
170pFront[1] *= -1
171frontR = occ.Box(tuple(pFront-fBox),
172 tuple(pFront+fBox))
173
174hAttach = carrierHeight*0.25*4
175wAttach = carrierWidth*0.8
176
177pAttach = np.array([carrierWidth*0.,carrierLength*0.,carrierHeight+hAttach*0.5])
178aBox = np.array([0.5*wAttach,0.5*wAttach,0.5*hAttach])
179
180attachment = occ.Box(tuple(pAttach-aBox),
181 tuple(pAttach+aBox))
182
183attachment.faces.Max((0,0,1)).name='attachment'
184
185#connection to sliders:
186sDim = 1.6*bDim
187sBoxes = []
188for ix in [-1,1]:
189 for iy in [-1,1]:
190 sBox = np.array([0.5*bDim,0.5*sDim,0.5*sDim])
191 psBox = np.array([ix*(carrierWidth*0.5-bDim*0.5),iy*(carrierLength*0.5-sDim*0.5),-sDim*0.5])
192 #print('psBox=',psBox)
193 box = occ.Box(tuple(psBox-sBox),
194 tuple(psBox+sBox))
195
196 cyl = occ.Cylinder(tuple(psBox-[0.5*bDim,0,0]), (1,0,0), rCable, bDim)
197 iName = 'joint_X'+str(ix)+'_Y'+str(iy)
198 cyl.faces[0].name = iName
199 interfaceNameListC += [iName]
200
201 sBoxes.append(box-cyl)
202
203
204geoC = occ.OCCGeometry(carrier+attachment+frontL+frontR
205 +sBoxes[0]+sBoxes[1]+sBoxes[2]+sBoxes[3])
206
207geoMeshC = geoC.GenerateMesh(maxh=maxhC,
208 curvaturesafety=curvaturesafety,
209 )
210
211meshC = Mesh(geoMeshC)
212
213femInterfaceC = FEMinterface()
214[bfM, bfK, fes] = femInterfaceC.ImportMeshFromNGsolve(meshC,
215 density=rho, youngsModulus=Emodulus, poissonsRatio=nu,
216 boundaryNamesList=interfaceNameListC,
217 meshOrder=2
218 )
219print('nNodesC=', femInterfaceC.NumberOfNodes())
220[boundaryNodesListC, boundaryWeightsListC] = femInterfaceC.GetBoundaryNodeSetsAsLists()
221femInterfaceC.ComputeHurtyCraigBamptonModes(boundaryNodesList=boundaryNodesListC,
222 nEigenModes=8,
223 excludeRigidBodyMotion=True,
224 boundaryNodesWeights=boundaryWeightsListC,
225 computationMode=HCBstaticModeSelection.RBE2)
226
227femInterfaceC.ComputePostProcessingModesNGsolve(fes, materialAlu)
228
229print('eigenfrequencies carrier (Hz):\n',np.round(femInterfaceC.GetEigenFrequenciesHz(),2),sep='')
230
231
232createFFRFObjectDictC = mbs.CreateFFRFReducedOrderObject(name='carrier',
233 referencePosition=[carrierWidth*0.5+bDim,
234 carrierLength*0.,
235 0.5*height+0.5*sDim],
236 #initialVelocity=[1,0,0],
237 femInterface=femInterfaceC,
238 stiffnessProportionalDamping=1e-4,
239 gravity=gravity)
240mCarrierAttachment = createFFRFObjectDictC['carrier:attachment']
241mCarrierSlidersLeft = [createFFRFObjectDictC['carrier:joint_X-1_Y-1'],
242 createFFRFObjectDictC['carrier:joint_X1_Y-1']]
243mCarrierSlidersRight = [createFFRFObjectDictC['carrier:joint_X-1_Y1'],
244 createFFRFObjectDictC['carrier:joint_X1_Y1']]
245
246
247if False: #activate to animate modes
248 mbs.Assemble()
249 from exudyn.interactive import AnimateModes
250
251 SC.visualizationSettings.nodes.show = False
252 #SC.visualizationSettings.view0.scene.showFaceEdges = True
253 SC.visualizationSettings.openGL.multiSampling=2
254 SC.visualizationSettings.openGL.lineWidth=2
255 SC.visualizationSettings.view0.window.renderWindowSize = [1600,1080]
256
257 objFFRFC = createFFRFObjectDictC['FFRFReducedOrderObjectDict']
258 nodeNumber = objFFRFC['nGenericODE2'] #this is the node with the generalized coordinates
259
260 SC.renderer.Start() #start graphics visualization
261 SC.renderer.SetModelView(zoom=2.011357,rotationVector=[-0.9999199,0.6224784,0.9588339],centerPoint=[0.7874073,1.334914,0])
262 AnimateModes(SC, mbs, nodeNumber, period=0.1, showTime=False,
263 renderWindowText='Eigenmodes visualization\n',
264 runOnStart=True)
265 import sys
266 sys.exit()
267
268
269#%%++++++++++++++++++++++++++++++++++++++++++++++++
270#add beams for rails:
271
272rhoA = 7800*rCable**2*np.pi
273
274fact = 1 #1e-2 #test with softer beam
275EA = rCable**2*np.pi*Emodulus*fact
276EI = rCable**4*np.pi/4*Emodulus*fact
277
278yBeam = (-0.5*width+thickness+1*bDim)
279
280ancfList = []
281for iy in [-1,1]:
282 p0 = np.array([0.5*bDimX,iy*yBeam,0.5*height])
283 p1 = p0 + [length-bDimX,0,0]
284
285 cable = ObjectANCFCable(physicsMassPerLength=rhoA,
286 physicsBendingStiffness = EI,
287 physicsBendingDamping = EI*0.001,
288 physicsAxialStiffness=EA,
289 physicsAxialDamping=EA*0.0001,
290 visualization=VObjectANCFCable(radius = rCable),
291 )
292
293 ancf=GenerateStraightLineANCFCable(mbs=mbs,
294 positionOfNode0=p0, positionOfNode1=p1,
295 numberOfElements=16, #converged to 4 digits
296 cableTemplate=cable, #this defines the beam element properties
297 massProportionalLoad = gravity,
298 #fixedConstraintsNode0 = [1,1,1, 0,1,1], #add constraints for pos and rot (r'_y,r'_z)
299 #fixedConstraintsNode1 = [1,1,1, 0,1,1], #add constraints for pos and rot (r'_y,r'_z)
300 )
301 #ancf=[cableNodeList, cableObjectList, loadList, cableNodePositionList, cableCoordinateConstraintList]
302
303 mBoxList = mFrameBox[0] if iy==-1 else mFrameBox[1]
304 nANCFnodes = len(ancf[0])
305 for i, marker in enumerate(mBoxList):
306 nANCF = ancf[0][i*(nANCFnodes-1)//2]
307 mANCF = mbs.AddMarker(MarkerNodePosition(nodeNumber=nANCF))
308 # mbs.CreateCartesianSpringDamper(bodyNumbers=[marker, nANCF],
309 # stiffness=[1e5]*3,
310 # damping=[2e3]*3)
311 mbs.AddObject(SphericalJoint(markerNumbers=[marker, mANCF],
312 visualization=VSphericalJoint(jointRadius=rCable*1.2)))
313
314 ancfList.append(ancf)
315
316
317
318#%%++++++++++++++++++++++++++++++++++++++++++++++++
319#sliding joint:
320sMarkerZ = 0.5*height+0.5*sDim-sDim*0.5 #ideal Z-position
321sMarkerY = (carrierLength*0.5-sDim*0.5)
322
323addSlidingJoint = True
324if addSlidingJoint:
325 for iANCF, markerList in enumerate([mCarrierSlidersRight,mCarrierSlidersLeft]):
326 ancf = ancfList[iANCF]
327 lElem = mbs.GetObject(ancf[1][0])['physicsLength']
328
329 for i, mSlider in enumerate(markerList):
330 pMarker = mbs.GetMarkerOutput(mSlider,exu.OutputVariableType.Position,
331 exu.ConfigurationType.Reference)
332 #print('pMarker=',pMarker)
333 if pMarker[1] < 0:
334 offsetY = -(sMarkerY - abs(pMarker[1]))
335 else:
336 offsetY = (sMarkerY - abs(pMarker[1]))
337 offsetZ = sMarkerZ- pMarker[2]
338 mbs.SetMarkerParameter(mSlider, 'offset', [0,offsetY,offsetZ])
339 pMarkerCorr = mbs.GetMarkerOutput(mSlider,exu.OutputVariableType.Position,
340 exu.ConfigurationType.Reference)
341 print('pMarkerCorr=',pMarkerCorr)
342
343 slidingCoordinateInit = pMarker[0] #X-coordinate
344 initialLocalMarker = int(slidingCoordinateInit/lElem) #second element
345
346
347 cableMarkerList = []#list of MarkerBodyBeamShape
348 offsetList = [] #list of offsets counted from first cable element; needed in sliding joint
349 offset = 0 #first cable element has offset 0
350 for item in ancf[1]: #create markers for cable elements
351 m = mbs.AddMarker(MarkerBodyBeamShape(bodyNumber = item))
352 cableMarkerList += [m]
353 offsetList += [offset]
354 offset += lElem
355
356 nodeDataSJ = mbs.AddNode(NodeGenericData(initialCoordinates=[initialLocalMarker,slidingCoordinateInit],numberOfDataCoordinates=2)) #initial index in cable list
357 slidingJoint = mbs.AddObject(ObjectJointSliding(markerNumbers=[mSlider,cableMarkerList[initialLocalMarker]],
358 constrainRotations=[0,0,0],
359 slidingMarkerNumbers=cableMarkerList, slidingMarkerOffsets=offsetList,
360 nodeNumber=nodeDataSJ))
361
362#++++++++++++++++++++++++++++++++++++++++++++++++
363#driving force and measurement:
364
365def UFspring(mbs, t, itemNumber, displacement, velocity, stiffness, damping, offset):
366 xDesired = SmoothStep(t, 0.2, 1.2, 0, length-carrierWidth-bDim*4)
367 k=1e5
368 d=5e3
369 forceX = (displacement[0]-xDesired)*k + velocity[0]*d
370 return [forceX,0,0]
371
372oGroundCSD = mbs.CreateGround(referencePosition=mbs.GetMarkerOutput(mCarrierAttachment,
373 exu.OutputVariableType.Position,
374 exu.ConfigurationType.Reference) )
375
376oCSD = mbs.CreateCartesianSpringDamper(bodyNumbers=[oGroundCSD, mCarrierAttachment],
377 springForceUserFunction=UFspring,
378 show=False)
379
380mbs.AddLoad(LoadForceVector(markerNumber=mCarrierAttachment, loadVector=[0,0,-2000*9.81]))
381
382sPos = mbs.AddSensor(SensorMarker(markerNumber=mCarrierAttachment, storeInternal=True,
383 outputVariableType=exu.OutputVariableType.Displacement))
384
385#++++++++++++++++++++++++++++++++++++++++++++++++
386mbs.Assemble()
387
388SC.visualizationSettings.nodes.show = False
389SC.visualizationSettings.markers.show = True
390SC.visualizationSettings.markers.drawSimplified = False
391SC.visualizationSettings.markers.defaultSize = 0.01
392
393if False: #activate to animate modes
394 mbs.Assemble()
395 from exudyn.interactive import AnimateModes
396
397 #SC.visualizationSettings.view0.scene.showFaceEdges = True
398 SC.visualizationSettings.openGL.multiSampling=2
399 SC.visualizationSettings.openGL.lineWidth=2
400 SC.visualizationSettings.view0.window.renderWindowSize = [1600,1080]
401
402 objFFRF = createFFRFObjectDict0['FFRFReducedOrderObjectDict']
403 nodeNumber = objFFRF['nGenericODE2'] #this is the node with the generalized coordinates
404
405 SC.renderer.Start() #start graphics visualization
406 SC.renderer.SetModelView(zoom=2.011357,rotationVector=[-0.9999199,0.6224784,0.9588339],centerPoint=[0.7874073,1.334914,0])
407 AnimateModes(SC, mbs, nodeNumber, period=0.1, showTime=False,
408 renderWindowText='Eigenmodes visualization\n',
409 runOnStart=True)
410 import sys
411 sys.exit()
412
413
414SC.visualizationSettings.contour.outputVariable = exu.OutputVariableType.StressLocal
415SC.visualizationSettings.contour.outputVariableComponent = -1
416
417SC.visualizationSettings.view0.window.renderWindowSize=[1200,800]
418SC.visualizationSettings.general.autoFitScene=False
419
420#SC.visualizationSettings.view0.scene.drawCoordinateSystem = False
421SC.visualizationSettings.openGL.lineWidth = 1
422SC.visualizationSettings.loads.show = False
423SC.visualizationSettings.openGL.light0.position=[2,2,10,1]
424
425SC.visualizationSettings.openGL.lightModelAmbient = [0.6,0.6,0.6,1]
426SC.visualizationSettings.openGL.multiSampling = 2
427SC.visualizationSettings.openGL.light0.shadow = 0.2
428SC.visualizationSettings.openGL.light1.shadow = 0.2
429SC.visualizationSettings.raytracer.numberOfThreads = 64
430
431SC.visualizationSettings.view0.camera.useRaytracer = False #set True for raytracing
432SC.visualizationSettings.raytracer.keepWindowActive= True
433SC.visualizationSettings.raytracer.advanced.searchTreeFactor = 8
434
435simulationSettings = exu.SimulationSettings()
436
437#simulationSettings.solutionSettings.writeSolutionToFile = False
438simulationSettings.solutionSettings.solutionWritePeriod = 0.02 #data not used
439simulationSettings.solutionSettings.sensorsWritePeriod = stepSize
440simulationSettings.timeIntegration.verboseMode = 1 #turn off, because of lots of output
441simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
442# simulationSettings.parallel.numberOfThreads = 4
443#simulationSettings.displayComputationTime = True
444#simulationSettings.displayStatistics = True
445
446simulationSettings.timeIntegration.numberOfSteps = int(endTime/stepSize)
447simulationSettings.timeIntegration.endTime = endTime
448simulationSettings.timeIntegration.newton.useModifiedNewton = True
449
450
451#visualize in Exudyn:
452SC.renderer.Start() #start graphics visualization
453
454SC.renderer.SetModelView(zoom=2.011357,rotationVector=[-0.9999199,0.6224784,0.9588339],centerPoint=[0.7874073,1.334914,0])
455SC.renderer.DoIdleTasks() #press space to continue
456mbs.SolveDynamic(simulationSettings)
457
458#SC.renderer.Stop() #safely close rendering window!
459
460mbs.PlotSensor(sPos, components=[0])
461mbs.PlotSensor(sPos, components=[2])
462
463mbs.SolutionViewer()