Dynamic simulation; Frequency domain response

[hamid401]

Hi all,
I'm trying to extract a system response in the frequency domain. I apply a step function of pressure on a system and do FFT on the velocity of a node. (instead on applying impulse load and doing FFT on displacement, I apply step load and do FFT on velocity). This should give me the magnitude and phase response of the system over a frequency range (depending on duration and sampling time).

I did similar simulations at Abaqus and Code-aster to compare:
In Code-aster, viscoelasticity is not defined in prony series for dynamic modelling, so I did dynamic linear elastic model with small damping (Rayleigh model) and harmonic simulations to confirm the dynamic response.
In Abaqus both hyperelastic and viscohyperelastic (VHE) models were modelled using the same coefficients as in FEBio.

Results:
(convergence tests were done to predict optimum frequency resolution (df) and frequency range (F) to obtain frequency response for F<1000 Hz)
Code-aster and Abaqus predict displacement magnitude very similar in both harmonic (steady state harmonic simulations) and dynamic with FFT on velocity.
FEBio does not predict any resonance at Frequencies <1000Hz for VHE model. Also for hyperelastic model, frequency response has quite noisy response and also does not show any anti-resonances as they are visible in Abaqus and Code-aster results.

Questions:
1. Is there any problem with my code in FEBio?
2. How much is dynamic modelling reliable in FEBio?
3. Is there any problem modelling viscoelastic materials in dynamic simulations?
4. Is there any plan to provide steady-state harmonic simulations for FEBio in future?
5. What is your suggestion to calculate frequency domain response of a system in FEBio?

Model description and graphs are in attached file (FFT_Report1 and 2) and FEBio file is also attached (I had to delete some lines of nodes, elements,BC and loading surface to reduce the size of the FEBio file)

[maas]

Hi Hamid,

Can you please send me the entire input file? If zipping it is still too big for upload, perhaps you can email it to me? steve dot maas at Utah dot edu.

Thanks,

Steve

[hamid401]

Hi Steve,
file is attached

[msbr]

Harmonic excitation, static simulation

Hi all,

any news on this topic?

We have a problem which is not too far away from the one described above, so I decided not to open a new thread but include it here:

In order to get the amplitudes of the particular solution of the displacement field of a harmonically excited solid at frequency omega, we substituted the original stiffness matrix K by an effective stiffness Keff=[K - (omega^2)M] in the source code, where M is the mass matrix. We then did a simple static simulation Keff*u0=F. By multiplying the amplitudes u0 with exp(j*omega*t) afterwards, we get the time-dependent solution again. Superposing the harmonic displacement field should therefore give a similar result to doing a full dynamic simulation.
The solid body consists of 10x10x1 uniformly sized Hex8 elements (edge length h=1).
The material is assumed to be linear-elastic, but non-homogeneous, as we ultimately want to detect stiff inclusions in the soft tissue by inverse analysis.
The nodes at the bottom are fixed in all directions, at the top the force acts in vertical direction.
We are using FEBio release 1.6.1 on Linux, with SuperLU or Skyline.

In general, the procedure works for very small omega, however as soon as we increase the frequency, the problem is diverging, giving a somehow strange error message of "negative jacobians". Of course, in order to obtain accurate results, the discretization has to be very fine, depending on the excitation frequency. But nevertheless, I don't really see why the Keff*u0=F cannot be solved for u0 already for frequencies that are only slightly above the smallest eigenvalues of M^(-1)K.
The divergence is neither affected by smaller time-step sizes, nor by smaller force amplitudes. Manipulating the convergence tolerances and max_refs, max_ups, max_retries etc. also only has minor influence. Besides, solving the full dynamic problem works also for much higher excitation frequencies.

My questions therefore are:
- How would manipulating the stiffness matrix have an impact on the element jacobians? I thought this was related to not suppressed rigid body modes or strange geometries of the element, or does this error message have a different meaning in the present case?
- Do you think that the direct solvers (following the LU decomposition approach) cannot handle the set of equations anymore when omega exceeds the smallest eigenvalues, i.e. when Keff is not positive definite anymore?
- If yes, has the Pardiso Solver any options to deal with matrices that are in general neither symmetric, nor positive definite? Or is it possible to include other iterative solvers in the FEBio-procedure?
- Or do you think that the problem is not related to the direct solver type at all? (Actually, Keff is invertible in MATLAB, but extracting the stiffness and importing in MATLAB is not practicable for finer discretizations)
- Apart from that, is there any possibility or workaround to do a modal analysis or eigenvalue analysis in FEBio?

I have included the input-file for the harmonic case (Keff u=F), note that we defined the force type as "force_harmonic" and included the excitation frequency omega in the load section. Apart from these two changes, the static case (Ku=F) is identical.
We changed only a small part in the source code, in the FESolidSolver.cpp file, function void "FESolidSolver::AssembleStiffness(vector<int>& en, vector<int>& elm, matrix& ke)": Before assembling m_pK into the global stiffness matrix, we constructed the element mass matrix, weighted the entries by the factor (omega^2) and subtracted it from the element stiffness matrix. As noted above, omega is read from the input-file.
In the present case, the simulation breaks for omega>5.5. I also attached K and Keff (as arrays) and a plot of the smallest eigenvalues of K, M^(-1)K and Keff (for omega=5.0). M was computed in MATLAB by (K - Keff) / (omega^2).

Hope that I explained everything that is relevant, and that you have some ideas why this is making so much trouble or even how we can solve/workaround that... Thanks a lot for help!

Best,
Chris

Short Update:
Solving Keff u=F in MATLAB works fine with LU decomposition for any arbitrarily high omega (using the matrices posted in the original post). The result was validated by comparison to a mode superposition of M^(-1).K, which gave the same result. Thus, I assume that the LU solver is just fine for that, but there is probably something "going wrong" after the assembly of the global stiffness matrix, meaning some sort of reformation/update/restructuring which makes the approach described above obsolete/useless.
What location in the code would you suggest for building Keff? I could also post the altered FESolidSolver.cpp if necessary.

Thanks for any response!

[kevinmattheusmoerman]

This approach seems interesting. Any news on this?

Kevin

[msbr]

Hi Kevin,

actually, we found that the problem is indeed not SuperLU, but the Quasi-Newton loop SuperLU is embedded in. So instead of calling SuperLU in each BFGS-iteration step for the calculation of displacement increments from the current residual and then updating the displacement solution with a line search algorithm, we now call SuperLU directly with the total displacement vector and the RHS-vector as arguments, which in our case mainly contains the concentrated nodal forces. In detail, the problem is the line search, because it searches for minima, and thus requires positive definite matrices. But in our case, where Keff=K-(w^2*M) is indefinite for large w, we are actually looking for saddle points. Solving the system of equations directly by means of the LU decomposition works for us.

At the moment we are investigating the accuracy of this approach for computing the amplitudes of the steady-state harmonic response, and the requirements concerning the mesh density. Btw, for these studies, we make use of your Gibbon toolbox, thanks a lot for providing it, works great!

Best,
Chris

[kevinmattheusmoerman]

Very interesting! I'd love to be able to simulate harmonic deformations with FEBio. We'd like to simulate magnetic resonance elastography, so like harmonic wave propagation in tissue. Do you think that is possible?

Let me know if there are any issues with GIBBON I can help with.


Kevin

[hamid401]

Any updates? Please let me know if I can do anything.
Hamid

[weiss]

Hi - A couple of comments:

Do you think that the direct solvers (following the LU decomposition approach) cannot handle the set of equations anymore when omega exceeds the smallest eigenvalues, i.e. when Keff is not positive definite anymore?

Yes, the stiffness matrix is likely becoming ill-conditioned, so the direct linear solve is returning a numerically inaccurate estimate for K^-1. Then the quasi-Newton method fails to converge.

actually, we found that the problem is indeed not SuperLU, but the Quasi-Newton loop SuperLU is embedded in. So instead of calling SuperLU in each BFGS-iteration step for the calculation of displacement increments from the current residual and then updating the displacement solution with a line search algorithm, we now call SuperLU directly with the total displacement vector and the RHS-vector as arguments, which in our case mainly contains the concentrated nodal forces. In detail, the problem is the line search, because it searches for minima, and thus requires positive definite matrices. But in our case, where Keff=K-(w^2*M) is indefinite for large w, we are actually looking for saddle points. Solving the system of equations directly by means of the LU decomposition works for us.

If you think the line search is causing the problem, you should be able to effectively disable it by adjusting the line search tolerance.

Most of the time I think these frequency domain analyses are done by finding the eigenvalues of the stiffness matrix. by subspace iteration. There are ways to extend this to handle higher modes more inexpensively. I don't think that we have this capability in FEBio right now, but it could be added.

Here is a relevant publication:

http://web.mit.edu/kjb/www/Principal_Publications/The_subspace_iteration_method_Revisited.pdf

The original publication was in 1973:

Bathe KJ, Wilson EL. Solution methods for eigenvalue problems in structural
mechanics. Int. J. Numer. Methods Eng. 1973;6:213?26.

Best regards,

Jeff

[weiss]

Hi - rather than replying to a 4 year old thread, please delete your response above and start a new thread with a new topic. This applies to the other old threads that you replied to as well