How to use Modal-superposition method for Harmonic analysis under base excitation

[Victor-Kao]

Hi everyone

Recently, I have been trying to perform the Modal-superposition method of Harmonic analysis under base excitation. The modal-superposition method (hope = MSUP) performs well if the excitation is force loading (F command); however, if the excitation is displacement (D command), the response will be 0 only. But if I use the (hropt = FULL) for base excitation, it can generate the correct response.
Therefore, I am wondering if there is any limitation of (hope = MSUP) that it can only work on force loading or if there is any method that I can use for base excitation

In the following, I wrote a simple SDOF system using COMBIN14 and MASS21 elements. And I tried three cases

1. First, using the force loading as excitation:

#from ansys.mapdl.core import launch_mapdl
import matplotlib.pyplot as plt
import numpy as np
import time 
#mapdl = launch_mapdl(nproc=4)

mapdl.clear()

## Prep7 #######################################################
mapdl.prep7()
mapdl.units("SI")  # SI - International system (m, kg, s, K).

# Create Node
mapdl.n(1,"0","0","0")
mapdl.n(2,"0","0","10")

# Define mass21 element 
mapdl.et(1, "MASS21", kop3=0)
mapdl.type(1)
mapdl.real(1)
mapdl.r(1,100,300,200)
mapdl.e(2)

mapdl.et(2, "COMBIN14", kop2=3)
mapdl.type(2)
mapdl.real(2)
k = (4*np.pi**2)*20000
mapdl.r(2,k,1000)
mapdl.e(1,2)

## SOLU #######################################################
mapdl.run("/SOLU")
mapdl.antype(2)
mapdl.modopt('LANB',3)
mapdl.mxpand('ALL','','','YES','','YES')
mapdl.outpr('ALL')
mapdl.outres("ALL","ALL")	
mapdl.outres("NSOL","1")
# Apply BC	
mapdl.d(1,"UZ",0)
mapdl.allsel()
mapdl.solve()
mapdl.finish()

mapdl.run("/SOLU")
mapdl.antype(3)
mapdl.hropt('MSUP')
mapdl.eqslv('SPARSE')
mapdl.harfrq("0.1","100")
mapdl.nsubst(int(100/0.2))
mapdl.dmprat("0.05")
# Apply BC	
mapdl.f(2,"FZ",1)
mapdl.d(1,"UZ",0)
mapdl.kbc(1)
mapdl.solve()
mapdl.finish()

mapdl.run("/SOLU")
mapdl.allsel()
mapdl.expass("ON")
mapdl.numexp("ALL")
mapdl.solve()
mapdl.finish()

## POST #######################################################
mapdl.post26()
a = mapdl.nsol("2",2,"U","Z","disp")
freq_values = np.unique(mapdl.post_processing.frequency_values)

# put the real and imaginary parts of acc_z_1F into unique mapdl array parameters
mapdl.vget("acc_z_1F_R", 2, '', 0)
mapdl.vget("acc_z_1F_IM", 2, '', 1)
real_part = mapdl.parameters["acc_z_1F_R"]
imag_part = mapdl.parameters["acc_z_1F_IM"]
#
plt.plot(freq_values, np.sqrt(np.power(real_part,2)+np.power(imag_part,2)), label="AMPL = sqrt(R^2+I^2)",linestyle="-.")
plt.legend()
plt.xlabel("freq (Hz)")
plt.ylabel("AMPL")
plt.show()

The plot looks like this, which is what I want:

2. Second, using the base excitation (D) with hropt = FULL, this can also generate the response I want :

#from ansys.mapdl.core import launch_mapdl
import matplotlib.pyplot as plt
import numpy as np
import time 
#mapdl = launch_mapdl(nproc=4)

mapdl.clear()

## Prep7 #######################################################
mapdl.prep7()
mapdl.units("SI")  # SI - International system (m, kg, s, K).

# Create Node
mapdl.n(1,"0","0","0")
mapdl.n(2,"0","0","10")

# Define mass21 element 
mapdl.et(1, "MASS21", kop3=0)
mapdl.type(1)
mapdl.real(1)
mapdl.r(1,100,300,200)
mapdl.e(2)

mapdl.et(2, "COMBIN14", kop2=3)
mapdl.type(2)
mapdl.real(2)
k = (4*np.pi**2)*20000
mapdl.r(2,k,1000)
mapdl.e(1,2)

## SOLU #######################################################
mapdl.run("/SOLU")
mapdl.antype(3)
mapdl.hropt('FULL')
mapdl.eqslv('SPARSE')
mapdl.harfrq("0.1","100")
mapdl.nsubst(int(100/0.2))
mapdl.dmprat("0.05")
# Apply BC	
mapdl.d(1,"UZ",1)
mapdl.kbc(1)
mapdl.solve()
mapdl.finish()


## POST #######################################################
mapdl.post26()
a = mapdl.nsol("2",2,"U","Z","disp")
freq_values = np.unique(mapdl.post_processing.frequency_values)

# put the real and imaginary parts of acc_z_1F into unique mapdl array parameters
mapdl.vget("acc_z_1F_R", 2, '', 0)
mapdl.vget("acc_z_1F_IM", 2, '', 1)
real_part = mapdl.parameters["acc_z_1F_R"]
imag_part = mapdl.parameters["acc_z_1F_IM"]
#
plt.plot(freq_values, np.sqrt(np.power(real_part,2)+np.power(imag_part,2)), label="AMPL = sqrt(R^2+I^2)",linestyle="-.")
plt.legend()
plt.xlabel("freq (Hz)")
plt.ylabel("AMPL")
plt.show()

The response looks like this:

3. Finally, I tried using the model-superposition method hropt = MSUP with the same BCs as the second case, but it gave me a 0 for the response;

#from ansys.mapdl.core import launch_mapdl
import matplotlib.pyplot as plt
import numpy as np
import time 
#mapdl = launch_mapdl(nproc=4)

mapdl.clear()

## Prep7 #######################################################
mapdl.prep7()
mapdl.units("SI")  # SI - International system (m, kg, s, K).

# Create Node
mapdl.n(1,"0","0","0")
mapdl.n(2,"0","0","10")

# Define mass21 element 
mapdl.et(1, "MASS21", kop3=0)
mapdl.type(1)
mapdl.real(1)
mapdl.r(1,100,300,200)
mapdl.e(2)

mapdl.et(2, "COMBIN14", kop2=3)
mapdl.type(2)
mapdl.real(2)
k = (4*np.pi**2)*20000
mapdl.r(2,k,1000)
mapdl.e(1,2)

## SOLU #######################################################
mapdl.run("/SOLU")
mapdl.antype(2)
mapdl.modopt('LANB',3)
mapdl.mxpand('ALL','','','YES','','YES')
mapdl.outpr('ALL')
mapdl.outres("ALL","ALL")	
mapdl.outres("NSOL","1")
# Apply BC	
#mapdl.d(1,"UZ",0)
mapdl.allsel()
mapdl.solve()
mapdl.finish()

mapdl.run("/SOLU")
mapdl.antype(3)
mapdl.hropt('MSUP')
mapdl.eqslv('SPARSE')
mapdl.harfrq("0.1","100")
mapdl.nsubst(int(100/0.2))
mapdl.dmprat("0.05")
# Apply BC	
mapdl.d(1,"UZ",1)
mapdl.kbc(1)
mapdl.solve()
mapdl.finish()

mapdl.run("/SOLU")
mapdl.allsel()
mapdl.expass("ON")
mapdl.numexp("ALL")
mapdl.solve()
mapdl.finish()

## POST #######################################################
mapdl.post26()
a = mapdl.nsol("2",2,"U","Z","disp")
freq_values = np.unique(mapdl.post_processing.frequency_values)

# put the real and imaginary parts of acc_z_1F into unique mapdl array parameters
mapdl.vget("acc_z_1F_R", 2, '', 0)
mapdl.vget("acc_z_1F_IM", 2, '', 1)
real_part = mapdl.parameters["acc_z_1F_R"]
imag_part = mapdl.parameters["acc_z_1F_IM"]
#
plt.plot(freq_values, np.sqrt(np.power(real_part,2)+np.power(imag_part,2)), label="AMPL = sqrt(R^2+I^2)",linestyle="-.")
plt.legend()
plt.xlabel("freq (Hz)")
plt.ylabel("AMPL")
plt.show()

It gave me the wrong result:

Please let me know if you have any suggestion.

Best,
Victor

[germa89]

Ok.. the following is comming from my old days at uni ... so double check this.

I believe it is a limitation of the MSUP itself, and not of the MAPDL MSUP implementation. You are using the code right, but for the wrong application.
I think that if you apply a displacement, and then you calculate the mode shapes (which is a needed step in MSUP to calculate the response as function of those mode shapes), the mode shapes where the displacement is applied gets cancelled. That's why when you apply a displacement in the base, the mode shapes do represent that node fully constrained (like a cantilever beam case for instance).

The approach on these cases is to introduce the "displacement" as an equivalent force by converting this displacement to an acceleration and then to force.
Check this article, it might give you some clues: https://simutechgroup.com/base-acceleration-in-harmonic-analysis-using-workbench/

[Victor-Kao]

Hi @germa89

Ok, thank you for your advice! I will do more research on this part.
Best,
Victor

[Victor-Kao]

Hi @germa89

One more question: I might also try Harmonic analysis using the Krylov method. Do you think this method can be used for such base excitation? or do I also need to consider transferring the displacement as an equivalent force?

Best,
Victor

[germa89]

@Victor-Kao

From the link you shared:

Similar to the full harmonic analysis, the frequency-sweep Krylov method uses full system matrices to compute the harmonic response.

Therefore, it is not using MSUP, hence it should be fine to use displacements.

[mikerife]

Hi @Victor-Kao & @germa89
MAPDL does not allow applied displacements with the D command with a mode supperposition harmonic analysis. The DVAL command is used, along with the MODCONT command. Also the D command is used in the modal solution to define the load identification number. The modal then harmonic solution would then be:

## SOLU #######################################################
mapdl.run("/SOLU")
mapdl.antype(2)
mapdl.modopt('LANB',4)
mapdl.mxpand('ALL','','','YES','','YES')
mapdl.outpr('ALL')
mapdl.outres("ALL","ALL")	
# mapdl.outres("NSOL","1")
# Apply BC
mapdl.modcont('', 'on')
mapdl.d(1,"UZ",1)
mapdl.allsel()
modal_solve = mapdl.solve()
mapdl.finish()

mapdl.run("/SOLU")
mapdl.antype(3)
mapdl.hropt('MSUP')
mapdl.eqslv('SPARSE')
mapdl.harfrq("0.1","100")
mapdl.nsubst(int(100/0.2))
mapdl.dmprat("0.05")
# Apply BC	
mapdl.dval(1, 'U', 1)
mapdl.kbc(1)
mapdl.solve()
mapdl.finish()

Mike

[Victor-Kao]

Hi @mikerife

Thank you for your suggestion! I tried your method; however, it has some error returns in the modal analysis.

from ansys.mapdl.core import launch_mapdl
import matplotlib.pyplot as plt
import numpy as np
import time 
mapdl = launch_mapdl(nproc=4)

mapdl.clear()

## Prep7 #######################################################
mapdl.prep7()
mapdl.units("SI")  # SI - International system (m, kg, s, K).

# Create Node
mapdl.n(1,"0","0","0")
mapdl.n(2,"0","0","10")

# Define mass21 element 
mapdl.et(1, "MASS21", kop3=0)
mapdl.type(1)
mapdl.real(1)
mapdl.r(1,100,300,200)
mapdl.e(2)

mapdl.et(2, "COMBIN14", kop2=3)
mapdl.type(2)
mapdl.real(2)
k = (4*np.pi**2)*20000
mapdl.r(2,k,1000)
mapdl.e(1,2)

## SOLU #######################################################
mapdl.run("/SOLU")
mapdl.antype(2)
mapdl.modopt('LANB',4)
mapdl.mxpand('ALL','','','YES','','YES')
mapdl.outpr('ALL')
mapdl.outres("ALL","ALL")	
# mapdl.outres("NSOL","1")
# Apply BC
mapdl.modcont('', 'on')
mapdl.d(1,"UZ",1)
mapdl.allsel()
modal_solve = mapdl.solve()
mapdl.finish()

mapdl.run("/SOLU")
mapdl.antype(3)
mapdl.hropt('MSUP')
mapdl.eqslv('SPARSE')
mapdl.harfrq("0.1","100")
mapdl.nsubst(int(100/0.2))
mapdl.dmprat("0.05")
# Apply BC	
mapdl.dval(1, 'U', 1)
mapdl.kbc(1)
mapdl.solve()
mapdl.finish()

## POST #######################################################
mapdl.post26()
a = mapdl.nsol("2",2,"U","Z","disp")
freq_values = np.unique(mapdl.post_processing.frequency_values)

# put the real and imaginary parts of acc_z_1F into unique mapdl array parameters
mapdl.vget("acc_z_1F_R", 2, '', 0)
mapdl.vget("acc_z_1F_IM", 2, '', 1)
real_part = mapdl.parameters["acc_z_1F_R"]
imag_part = mapdl.parameters["acc_z_1F_IM"]
#
plt.plot(freq_values, np.sqrt(np.power(real_part,2)+np.power(imag_part,2)), label="AMPL = sqrt(R^2+I^2)",linestyle="-.")
plt.legend()
plt.xlabel("freq (Hz)")
plt.ylabel("AMPL")
plt.show() 

It returns the error:

I assumed this error came from a lack of enough BCs since it mentioned Please check for an insufficiently constrained
model., Therefore, I add the BCs at the UX and UY dirs of node 2 ( I used the COBIN14 element with SDOF here; I assumed I didn't need to do this ).

## SOLU #######################################################
mapdl.run("/SOLU")
mapdl.antype(2)
mapdl.modopt('LANB',4)
mapdl.mxpand('ALL','','','YES','','YES')
mapdl.outpr('ALL')
mapdl.outres("ALL","ALL")	
# mapdl.outres("NSOL","1")
# Apply BC
mapdl.modcont('', 'on')
mapdl.d(1,"UZ",1)
mapdl.d(2,"UX",0)
mapdl.d(2,"UY",0)
mapdl.allsel()
modal_solve = mapdl.solve()
mapdl.finish()

mapdl.run("/SOLU")
mapdl.antype(3)
mapdl.hropt('MSUP')
mapdl.eqslv('SPARSE')
mapdl.harfrq("0.1","100")
mapdl.nsubst(int(100/0.2))
mapdl.dmprat("0.05")
# Apply BC	
mapdl.dval(1, 'U', 1)
mapdl.kbc(1)
mapdl.solve()
mapdl.finish()

It can work; however, the result is not what I expected as an FRF of the SDOF system:

Currently, I am using the 2021 R2 version of the ANSYS; I am not sure whether the version difference matters or not...
Please let me know if you have any suggestions.

Best,
Victor

[mikerife]

Hi @Victor-Kao turn off pivot checking prior to the modal solve with

mapdl.pivcheck('off')

[Victor-Kao]

Hi @mikerife

Thank you very much for the helping! The code can work now, however, the one last problem is that the amplitude of the FRF looks incorrect. Attached is the comparison of using "FULL" and "MSUP.", as you can see, it seems like the damping ratio didn't apply correctly for the MSUP method (I used the mapdl.dmprat( ) to apply damping ratio)

Do you know what the problem might be with this issue?

Here is the code I used for MSUP:

#from ansys.mapdl.core import launch_mapdl
import matplotlib.pyplot as plt
import numpy as np
import time 
#mapdl = launch_mapdl(nproc=4)

mapdl.clear()

## Prep7 #######################################################
mapdl.prep7()
mapdl.units("SI")  # SI - International system (m, kg, s, K).

# Create Node
mapdl.n(1,"0","0","0")
mapdl.n(2,"0","0","10")

# Define mass21 element 
mapdl.et(1, "MASS21", kop3=0)
mapdl.type(1)
mapdl.real(1)
mapdl.r(1,100,300,200)
mapdl.e(2)

mapdl.et(2, "COMBIN14", kop2=3)
mapdl.type(2)
mapdl.real(2)
k = (4*np.pi**2)*20000
mapdl.r(2,k,1000)
mapdl.e(1,2)

## SOLU #######################################################
mapdl.run("/SOLU")
mapdl.antype(2)
mapdl.modopt('LANB',4)
mapdl.mxpand('ALL','','','YES','','YES')
mapdl.outpr('ALL')
mapdl.outres("ALL","ALL")	
# mapdl.outres("NSOL","1")
# Apply BC
mapdl.modcont('', 'on')
mapdl.d(1,"UZ",1)
mapdl.dmprat(0.05)
mapdl.allsel()
mapdl.pivcheck('off')
modal_solve = mapdl.solve()
mapdl.finish()

mapdl.run("/SOLU")
mapdl.antype(3)
mapdl.hropt('MSUP')
mapdl.eqslv('SPARSE')
mapdl.harfrq("0.1","100")
mapdl.nsubst(int(100/0.2))
mapdl.dmprat(0.05)
# Apply BC	
mapdl.dval(1, 'U', 1)
mapdl.kbc(1)
mapdl.solve()
mapdl.finish()

mapdl.run("/SOLU")
mapdl.allsel()
mapdl.expass("ON")
mapdl.numexp("ALL")
mapdl.solve()
mapdl.finish()
## POST #######################################################
mapdl.post26()
a = mapdl.nsol("2",2,"U","Z","disp")
freq_values = np.unique(mapdl.post_processing.frequency_values)

# put the real and imaginary parts of acc_z_1F into unique mapdl array parameters
mapdl.vget("acc_z_1F_R", 2, '', 0)
mapdl.vget("acc_z_1F_IM", 2, '', 1)
real_part = mapdl.parameters["acc_z_1F_R"]
imag_part = mapdl.parameters["acc_z_1F_IM"]
#
plt.plot(freq_values, np.sqrt(np.power(real_part,2)+np.power(imag_part,2)), label="AMPL = sqrt(R^2+I^2)",linestyle="-.")
plt.legend()
plt.xlabel("freq (Hz)")
plt.ylabel("AMPL")
plt.show()

And the solver for the "FULL" harmonic analysis:

## SOLU #######################################################
mapdl.run("/SOLU")
mapdl.antype(3)
mapdl.hropt('FULL')
mapdl.eqslv('SPARSE')
mapdl.harfrq("0.1","100")
mapdl.nsubst(int(100/0.2))
mapdl.dmprat(0.05)
# Apply BC	
mapdl.d(1,"UZ",1)
mapdl.kbc(1)
mapdl.solve()
mapdl.finish()

Thank you again!
Best,

Victor

[mikerife]

@Victor-Kao Use the QR Damped modal solver. See section 4.7.1 here for more information:
https://ansyshelp.ansys.com/public/account/secured?returnurl=/Views/Secured/corp/v242/en/ans_str/Hlp_G_STR4_MODESUPER.html%23aQcr1c5lcd