Custom field output

[Shah]

Hi all,

I am planning to generate a derived field output for Von-Mises stress response due to random vibration from harmonic results. To do that, I will need to process multiple frequency stress outputs, which I do not think can be done using the GUI at this point.

So I am thinking of accomplishing it externally to GUI, using a Python script by following the steps below:

• Read in the .frd file and process the multi-frequency element stress tensors (using both real and imaginary results)
• Create a new fictitious output step and add the new results (scalar Von-Mises values) as a standard output (e.g.: STRESS SXX).
• Open the processed frd file in PrepoMax to visualize the response.

I am wondering if anyone has the experience in such custom postprocessing and can advice if there’s a better way to achieve this before I go too deep.

[FEAnalyst]

How about this approach: Random vibration analysis of a cantilever beam - #2 by FEAnalyst ?

Might be good to ask on the CalculiX forum too. Some guys there dealt with similar issues. Here’s an example (response spectrum analyses): Shock Response Spectrum - Is there a workaround? - CalculiX (official versions are on www.calculix.de, the official GitHub repository is at https://github.com/Dhondtguido/CalculiX).

Also, it seems that Code_Aster supports random vibration analyses.

[Shah]

Hi Jakub,

Thanks a lot for pointing me to the reference random vibration problems with base excitation. I noticed the Abaqus example of a cantilever beam that you had mentioned in a related post. And I gave it a quick try to see if the FRF (freq response function) method I was planning would generate comparable results.

So I used the steady state dynamic analysis to generate FRF with the unit base motion, and used it to convert the input PSD to response acceleration RPSD using the linear transformation:

RPSD_accel = PSD*g^2*FRF^2

Then, I converted it to the displacement RPSD using:

RPSD_disp = RSPD_accel / w^4

where w = angular frequency in rad/s. I multiplied the RPSD_disp by 3 to get the 3 sigma response (to compare with the reference solution).

The following is the plot of the FRF curve for first two modes:

And here’s the comparison of calculiX to the reference solution displacement PSD:

As you can see, the peaks and the levels in between the two adjacent modes are very similar. But interestingly, the calculiX based response is overpredicting the RPSD levels relative to the reference at the beginning and at the end of the spectrum before the first mode and after the second mode. I wonder where this difference is coming from. From what I understand, the FRF is expected to be close to unity at the low freqs, and based on that the equations above are expected give such high levels near f=20Hz. Is there a special treatment applied to the ends of the curve in random vibration analysis that I may be missing? Any insights on this?

[FEAnalyst]

@Shah

So I used the steady state dynamic analysis to generate FRF with the unit base motion

Can you share the input file ? I would check the result in Abaqus.

@Shah

Is there a special treatment applied to the ends of the curve in random vibration analysis that I may be missing?

Mainly the concept of setting the number of data points before and after (and between) the eigenfrequencies but this is also used in CalculiX. Also, in Abaqus it’s recommended to use residual modes to avoid modal truncation errors but it’s not done by default and not available in CalculiX.

[Shah]

Sure, here’s the link to the input file, prepared using Prepomax. I applied unit base displacement motion using *boundary, and obtained the displacement response at the tip (which is my FRF). The modal damping is 0.1. The steady state dynamics analysis uses the first 10 eigenmodes.

Dropbox

To get such a low response at freq=20Hz in the order of RPSD = 7e-11 from abaqus, I would need to set the FRF to zero at that location. So it looks like the FRF is being tapered off in some ways at the beginning and at the higher end.

[FEAnalyst]

And then how do you postprocess it to obtain FRF ? I get this from your file:

Also, I wanted to submit it in Abaqus but it doesn’t allow non-zero BCs in mode-based SSD so it’s necessary to either use *BASE MOTION (CalculiX also has it, I thought you were using it instead of *BOUNDARY) or direct SSD procedure but then modal damping is not supported and other forms of damping have to be used instead.

[Shah]

I am using the magnitude of the tip displacement as FRF, I think you plotted the real displacement outputs. The magnitude is basically giving the max displacement amplitude at the given frequency.

In Abaqus, base motion with displacement input might be OK if you want to try it. I have tested both boundary and base motions in calculix with the same FRF outcome for the tip.

[FEAnalyst]

Ok, I see. This is the FRF that I got from Abaqus analysis with base motion:

FRF abq.zip (12.2 KB)

[Shah]

Thanks a lot Jakub, this sheds a lot more light to what may be going on here. Just so I understand, the abaqus FRF is also the tip displacement magnitude right (the same way calculix output is generated)?
Here’s the comparison between the two. It is interesting to see these curves together, the fact that the FRF level approaches zero at low and high freqs, and relatively higher level in abaqus curve between the two modes. I looks like Abaqus is generating the outputs differently than calculix. I need to get to the bottom of this, because this seems to be the source of the difference.

Using this abaqus curve, I got reasonably similar levels to the reference w/o amplifying the results by 3 (it looks like the output is not 3sigma response after all). So, I removed the factor 3 from calculix response as well in the plot below. It looks like the approach above to predict the RPSD is good, but now I need to figure out where the difference in steady state dynamics analysis is coming from.

[FEAnalyst]

@Shah

Just so I understand, the abaqus FRF is also the tip displacement magnitude right (the same way calculix output is generated)?

Yes, I used the same input file (just adapted for Abaqus) and selected magnitude in complex form settings. The only significant difference is the use of base motion:

*AMPLITUDE, NAME=amp1
10,1
1000,1
…
*BASE MOTION, DOF=2, AMPLITUDE=amp1, TYPE=DISPLACEMENT

[Shah]

Thank you for the clarification. I think I got to the bottom of this. It looks like abaqus prints out the relative displacement while calculix outputs the global displacement. To convert calculix displacement outputs to relative displacement, I would need to subtract input displacement from the tip displacement in complex form to account for the phase difference and then determine the magnitude. For example, the output is given in complex form as:

uy = uy_{real} + i*uy_{imag}

And the input is:

uy0 = uy0_{real} + i*uy0_{imag}

where uy0_{real} =1 and uy0_{imag}=0. Then the relative displacement would be:

uy_{relative} = uy_{real} - 1 + i*uy_{imag}

The magnitude of this would give the same curve as abaqus FRF curve as below:

By transforming the input PSD to output RPSD using this FRF, we can get the same response as the abaqus, both of which are comparable to the reference:

Thanks a lot once again for your help to resolve the difference. I feel more confident to apply this FRF technique to determine the von Mises response due to random vibration given the input PSD.

[FEAnalyst]

@Shah

It looks like abaqus prints out the relative displacement while calculix outputs the global displacement.

Yes, that’s right. From Abaqus documentation about mode-based SSD:

The standard output variables U, V, A, and the variable PU listed above correspond to motions relative to the motion of the primary base in a mode-based analysis. Total values, which include the motion of the primary base, are also available:

TU - Magnitude of all components of total displacement/rotation at a node.
…

I’m glad you’ve managed to solve it.

[lucas_bueno]

I have little to no experience in FRF analysis, but it is always nice to see that CCX can handle this type of analysis as well. Nice comparison

[Shah]

Hi, I finally pulled together a script to estimate the von-Mises RMS output and to map the results to the mesh via the frd file. It uses the same FRF approach above but for the von-Mises stress, although I may be over-simplifying the problem without seeing the whole picture. I mapped the stresses to the energy output in the frd file as it appeared to be the simplest scalar output type, but as expected, the unit is incorrect, the results are supposed to be in [MPa].

This type of output is not provided in the Abaqus benchmark study at the link below, but it would be very important to know RMS stresses for fatigue life evaluations.

https://classes.engineering.wustl.edu/2009/spring/mase5513/abaqus/docs/v6.6/books/bmk/default.htm?startat=ch01s04ach47.html

Jakub, if you still have the Abaqus input for this 3D model, it would be nice if you could drop such von-Mises RMS contour plot here for a qualitative comparison. It would be interesting to see if this contour plot is in the right ballpark.

[FEAnalyst]

@Shah

Jakub, if you still have the Abaqus input for this 3D model, it would be nice if you could drop such von-Mises RMS contour plot here for a qualitative comparison.

Sure, here’s a contour plot of RMS Mises output from random response analysis:

[Shah]

Thanks a lot Jakub, it is great to see the good agreement between the two. The FRF approach seems to be working just fine for the RMS Mises stresses too. This should be because the Mises stress output is a linear function of the input base acceleration in harmonic analysis. I noticed Miles’ approach also does a very good job capturing the RMS stresses for this type of models when the major mode is the flex mode that takes the similar shape to the displacement response under the static gravity load.

Thanks again for your help with the verification.

[dirkmunro]

Hi Shah

I am interested in this method. (Also in an attempt to make sense of Ansys results.)

Can you perhaps ellaborate and/or share input files? I just cant seem to reproduce the RMS Mises stresses exactly in the example from above.

Thank you very much in advance.

**

As I wrote the message, I think I may have figured it out:

The Mises stress which is integrated over the frequency domain, has to be computed from the magnitude of both the real and the imaginary stress quantities (not just the real part).

[Shah]

Hi Dirkmunro,

It looks like you already figured it out. In case you want to compare with my approach, here are the main python code snippets I used.

For each element, tensor and tensori are the stress tensors for the real and imaginary component of the stresses. The vonMises stress is vmRow (element level).

  tensor = np.array([[strsRow.SXX, strsRow.SXY, strsRow.SZX],
                     [strsRow.SXY, strsRow.SYY, strsRow.SYZ],
                     [strsRow.SZX, strsRow.SYZ, strsRow.SZZ]])

  tensori = np.array([[strsRowi.SXX, strsRowi.SXY, strsRowi.SZX],
                      [strsRowi.SXY, strsRowi.SYY, strsRowi.SYZ],
                      [strsRowi.SZX, strsRowi.SYZ, strsRowi.SZZ]])

  strsCmplx = tensor + 1j * tensori
  S1, S2, S3 = np.linalg.eigvals(strsCmplx)
  tau = np.array([S1 - S2, S2 - S3, S3 - S1])
  vmRow = sqrt(0.5) * np.linalg.norm(tau)

Then, in another routine, this vmRow comes in as eldat.Svm, which is a vector that carries VM stress for all frequencies. A few lines are to convert the FRF per acceleration to per g (gravity). Then, as you mentioned, the integration over all frequencies is carried out to covert the entire spectrum to a single random VM stress response.

    eldat["FRFstrsPerDisp"] = eldat.Svm  
    eldat["w"] = eldat.freq*2*pi   # angular freq, rad/s
    eldat["FRFstrsPeraccel"] =  eldat["FRFstrsPerDisp"]/eldat.w**2.
    eldat["FRFstrsPerG"] = eldat["FRFstrsPeraccel"]*9810.
    eldat["ASD"] = ASD(eldat.freq)
    eldat["Gstrs_out"] = eldat["FRFstrsPerG"]**2*eldat["ASD"]
    sum = np.trapezoid(eldat.Gstrs_out, x=eldat.freq)
    vmStrsRandom = sqrt(sum)

[dirkmunro]

Indeed, but thank you very much. I apply a base motion (typpe=acceleration) as 1 G (9810 mm/s^2), and take the square of the stress field at a particular frequency, to get the ‘transfer function‘ in MPa^2 / G^2. This gets multiplied the PSD value (G^2 / Hz ) at that frequency, which is a number in MPa^2/ Hz, which, when integrated over frequencies, gives MPa^2.

Note that there exists also the so-called Segalman-Reese algorithm (link below). With this, the transfer function is analytical, and the Mises RMS is calculated based only on the eigenmodes and frequencies. (Not a subsequent harmonic step from which the transfer function is obtained.) But I have not been able to make this work… Step 1 would be to reproduce the results of the harmonic step, with only the eigenmodes and the analyttical transfer function.

https://www.osti.gov/servlets/purl/753406