Trouble during a base excitation harmonic analysis of a cantilever beam

[JoiKjerulf]

Hey all, I'm trying to write a script to automate the analysis of the frequency response of cantilever beams to a base excitation and am running into issues with the data that I get as it does not match the manual analysis I've done in Workbench. I'm doing this in a jupyter workbook, with MAPDL version 25.1 and Ansys student version 2025 R1.

I'm using the mode-superposition method as described here, and so I am using the QRDAMP modal analyser, and DVAL to apply the displacement in the harmonic analysis step. I have looked at and tried to see what I'm doing wrong from question 3508 but I can't see it.

The beam setup and meshing is done as follows:

from ansys.mapdl.core import launch_mapdl
import numpy as np
import matplotlib.pyplot as plt
from ansys.math.core.math import AnsMath

mapdl = launch_mapdl(set_no_abort=True)

mapdl.clear()

mapdl.run("/FILNAM,MYJOB")           # change MYJOB to your desired job name
mapdl.run("/TITLE,Enforced Motion Example")
mapdl.prep7()



# All units in (m, Kg, s)
LENGTH = 6
RADIUS = 0.3
RHO = 7850
E = 2e11
NU = 0.27

mapdl.clear()
mapdl.prep7()
mapdl.cyl4(0,0,RADIUS,'','','',LENGTH)

#mapdl.vplot('all')
mapdl.lesize("ALL", 0.15, layer1=1)

mapdl.mp('ex',1,E)
mapdl.mp('nuxy',1,NU)
mapdl.mp("DENS" , 1, RHO)

mapdl.et(1,'SOLID186')
mapdl.mshape(1, "3D")
mapdl.mshkey(0)
mapdl.vmesh('all')
#mapdl.eplot()

Then I do the modal analysis like so

mapdl.nsel("S", "LOC", "Z", 0.0)     # e.g. select nodes at X=0
mapdl.cm("DRIVE_END", "NODE")            # component BASE1 = those nodes

mapdl.finish()

mapdl.run("/SOLU")
mapdl.antype("MODAL")
mapdl.modopt("QRDAMP", 10,"","","OFF")             # Block Lanczos, 10 modes
mapdl.mxpand("ALL","","","YES","","YES")
mapdl.outpr('ALL')
mapdl.outres("ALL","ALL")	

# Activate enforced-motion mode calculation
mapdl.modcont('', 'on')
mapdl.pivcheck('off')

# Apply enforced‐motion constraints
mapdl.cmsel("S", "DRIVE_END", "NODE")    # select component BASE1
mapdl.d("ALL", "ALL", 0)
mapdl.d(node="ALL", lab="UZ", value=1,)

mapdl.allsel("ALL")
mapdl.solve()

Which does give me correct resonant frequencies.

1 11.74488427782
2 11.74631933448
3 71.29438343423
4 71.30377615882
5 131.9752568233
6 190.6455019892
7 190.6704809331
8 210.7601646460
9 352.2418386847
10 352.2879443978

Then I do the harmonic analysis like so:

mapdl.finish()

mapdl.run("/SOLU")               
mapdl.antype("HARMIC")             
mapdl.hropt("MSUP", 10)        
mapdl.run("HROUT") #,CLUSTER,ON,4")    
mapdl.run("LVSCALE,1,1.")          
#mapdl.f("...", "FY", "...")    
mapdl.harfrq(1, 500)          
mapdl.dmprat(0.02)          
#mapdl.mdamp(1, 0.02)          
mapdl.nsubst(250)              
mapdl.cmsel("S", "DRIVE_END", "NODE")
mapdl.dval(1, 'U', 1)
mapdl.kbc(1)                        
#mapdl.save()                    
mapdl.allsel()
mapdl.solve()       

Then for the post processing, which I believe is my issue, although I haven't been able to verify if it's the harmonic analysis which is the problem, or if it's the post processing.

mapdl.finish()

mapdl.post26()
mapdl.file(mapdl.jobname, "rfrq")           # FILE,Jobname,RFRQ

mapdl.nsel("S", "LOC", "Z", LENGTH)
mapdl.get("TIPNODE", "NODE", 0, "NUM", "MIN")


mapdl.get("NCNT", "NODE", 0, "COUNT") #This was just to verify I'm selecting any nodes

tip_node = int(mapdl.parameters["TIPNODE"])
print(tip_node)
uz_cx = mapdl.nsol(2, tip_node, "U", "Z")

freq = mapdl.post_processing.frequency_values  

mapdl.finish()

uz_cx = uz_cx.flatten()        
npts   = uz_cx.size            
fmin, fmax = 1.0, 500.0

freq = fmin + np.arange(fmin, fmax, (fmax-fmin)/250)
amp  = np.abs(uz_cx)

This results in a strange graph. I'd have expected peaks around the resonance frequencies and not this.

I'm really hoping this is just a case of me misunderstanding something trivial, thanks for reading if you made it this far :-)

[mikerife]

Hi @JoiKjerulf took me a while to figure this out - I had a typo and so was still getting the same result you were and it took me forever to realize what I mistyped! Any way, uz_cx is a complex valued MAPDL array. To get the amplitude try the following instead (I changed the jobname to 'file'):

mapdl.post26()
mapdl.file('file', "rfrq")

mapdl.nsel("S", "LOC", "Z", LENGTH)
mapdl.get("TIPNODE", "NODE", 0, "NUM", "MIN")
mapdl.nsol(2, int(mapdl.parameters["TIPNODE"]), "U", "Z")

# get the real and imaginary values of variable 2 into separate mapdl arrays
mapdl.vget('uz_R', 2, '', 0)
mapdl.vget('uz_I', 2, '', 1)
# calculate amplitude
amp  = np.sqrt(mapdl.parameters["uz_R"]**2 + mapdl.parameters["uz_I"]**2)

fmin, fmax = 1.0, 500.0
freq = fmin + np.arange(fmin, fmax, (fmax-fmin)/250)

plt.plot(freq, amp)
plt.grid(True)
plt.show()

Then the plot will be as expected:

Mike

[germa89]

btw @mikerife , I would have expected the following code to plot the same as your graph without having to square anything:

mapdl.post26()
mapdl.file('file', "rfrq")

mapdl.prcplx(1) # Changing to Amplitude-phase

tip_node_amp = mapdl.nsol(9, tip_node, "U", "Z")

plt.plot(freq, tip_node_amp)
plt.grid(True)
plt.show()

[germa89]

Ok.. it seems that PRCPLX only works with PRVAR. But we do not use that to extract the variable values. We converted the stored variable (always in real-imaginary numbers) to an array (directly) using VGET, hence this is not "listing" per se.
The conversion should be done behind the scenes.

[germa89]

Issue created (another one): #3932

[JoiKjerulf]

Thank you for taking a look at this for me, this does indeed work and was just me misunderstanding what I was working with! You guys are great :-)

[germa89]

Kind remainder... you can do the following in one line:

mapdl.get("TIPNODE", "NODE", 0, "NUM", "MIN")
tip_node = int(mapdl.parameters["TIPNODE"])

with

>>> tip_node  = mapdl.get("TIPNODE", "NODE", 0, "NUM", "MIN")
>>> print(tip_node)
94

Or

tip_node  = mapdl.get_value("TIPNODE", "NODE", 0, "NUM", "MIN")

Same with vget which now returns automatically the array:

>>> uz_R = mapdl.vget('uz_R', 2, '', 0)
>>> print(uz_R)
array([[  1.00025719],
       [  1.00071431],
       [  1.0014005 ],
       [  1.0023166 ],
       [  1.00346446],
       ...
       ])

The command mapdl.nsol also returns an array:

>>> tip_node_nsol = mapdl.nsol(2, tip_node, "U", "Z")
>>> print(tip_node_nsol)
[[  1.00025719]
 [  1.00071431]
 [  1.0014005 ]
 [  1.0023166 ]
 [  1.00346446]
 [  1.0048403 ]
 ...
 ])

but only the amplitud or real part of the given variable, since mapdl.nsol does not pass any kcplx value to VGET.

>>> tip_node_nsol = mapdl.nsol(2, tip_node, "U", "Z")
>>> uz_R = mapdl.vget('uz_R', 2, '', 0)
>>> np.allclose(uz_R, tip_node_nsol) 
True

[JoiKjerulf]

That's very useful, thank you!