BeamDyn - Solution does not converge after the maximum number of iterations

[JelmerPolman]

Dear all,

I have some problems with my BeamDyn model. The following error keeps popping up at some simulations: "FAST_Solution:FAST_AdvanceStates:B1:BD_GA2:BD_DynamicSolutionGA2:Solution does not converge after the maximum number of iterations". I have reduced the time step to 0.0001 seconds, which results in unreasonable simulation times.

In order to troubleshoot I simulated the following testcase. Blade 1 is pitched 90 degrees and is in the three o'clock position and loaded with 1 g gravity. The rotor does not turn and AeroDyn is turned of so that only the deflection due to gravity is simulated. After a few time steps the simulation aborts.

QuasiStaticInit is set to true. One can see that the blade deflections are indeed static (1.3 m), but the root bending moment in x-direction (although very small at 17 N.m) is not.

The only thing that really helps to prevent this error is reducing the element order of BeamDyn. An element order of 5 for example runs without a problem. However, this order does not provide accurate results. I did a convergence study and calculated a reference solution with Ansys. From a comparison it follows that I need a element of order of at least 11 to get correct results.

If we look at the blade deflections, it can be clearly be seen that a low element order does not give accurate results. However, a high element order leads to a Runge's like oscillation in the torsion of the blade tip.

What can be the cause for OpenFAST to abort, althoug almost al channels are static? What steps can I take to reduce the oscillations in the blade torsion. Could this also be the cause for the needed very small time step?

I have attatched the input files here:

Simulation_files.zip

Best regards, Jelmer Polman

[andrew-platt]

Hi Jelmer,

Thanks for attaching the files. I have been working to remedy a couple of known issues that may be contributing to what you are seeing here. I plan on issuing a pull request to within a few weeks that will include some of the following fixes

1. Improved internal loads calculations (the methodology in OpenFAST 2.2.0 is incorrect)
2. Improve the calculations for reference curvature (this is causing incorrect internal strains at any orientation other than reference orientation)
3. Update to the cubic spline algorithm (this has been leading to incorrectly capturing the curvature and twist at the root and tip of the blade)
4. Improved reporting to the user on how the shape functions and cubic spline match the given keypoint line

Issues 2 and 3 are known to lead to non-convergence. If you would like to test the updated code that fixes these issues, you can access it on my branch at https://github.com/andrew-platt/openfast/tree/bug/BD_LoadsOutput (these updates are not ready for merging into the main code yet as we are still testing).

Regards, Andy

[JelmerPolman]

Hi Andy,

Thank you for sharing your code.

I have compiled it and now I am running a convergence study and some reference testcases. This will take some time and I will update you on the results as soon as possible!

Best regards, Jelmer Polman

[JelmerPolman]

Hi Andy,

wtih your version, the QuasiStatic initialization fails on the dynamic non-convergence. If I disable it and run a time simulation with a heavily damped blade, I was still able to do the convergence study.

I have found that it doesn't converge to the Ansys reference solution. From order 26 and onwards, the simulation even fails on the dynamic non-convergence.

The cross-sectional matrices at the quadrature points remained the same. However, the quadrature points changed, leading to different full stiffness matrix and mass matrix and thus different results. I thought that the quadrature points are set equal to the keypoints if the Trapezoidal method is chosen. How are the quadrature points defined in your version?

I add the BeamDyn summary files for both versions:

Test_07_RefBlade_OpenFAST_v2_2_0_Nel1_P11.BD1.txt

Test_07_RefBlade_OpenFAST_vandy_Nel1_P11.BD1.txt

Best regards, Jelmer Polman

[andrew-platt]

Hi Jelmer,

Thanks for your update and comparison. I agree that there is something obviously wrong with the latest update on my branch. I have not had a chance to look into this in detail, but hope to in the next couple of weeks.

I'm not surprised that order_elem of 26 or higher does not converge. There simply aren't a sufficient number of quadrature points in the mass and stiffness matrices to fully resolve the finite element points (see #183).

There might be some confusion here regarding the keypoints. They keypoints are used to define the curvature line of the blade. The distance along this curved line is then supposed to be used used for setting the location of all the mass and stiffness matrices (the blade definition file sets those according the their eta value). So it is possible to define the keypoint line with as few as 4 points, but still have a very large number of quadrature points.

The problem we discovered is that the quadrature points were being set by the z-location relative to the tip location instead of at the eta location along the blade curvature. This is incorrectly defining it for the geometrically exact beam theory (GEBT). For straight beams, it isn't an issue since the the curve defining the blade is a straight line along the z-axis. So, the latest set of commits was an attempt to correct this, but apparently something is now incorrect.

I really appreciate your feedback on this. Your comparison with the ANSYS results is extremely valuable in confirming that something is incorrect.

Regards, Andy

[JelmerPolman]

Hi Andy,

thank you for your update.

If you want, I would be more than willing to check your code against ANSYS. Just let me now if you commit a possible fix.

Regards, Jelmer

[JelmerPolman]

Hi Andy,

the research I am doing is related to a code-to-code comparison with other wind turbine simulation tools.

My BeamDyn model compares good with the results of other tools for the flapwise loading testcase as described in my first post, apart from a small difference in the torsional rotation at the blade tip:

However, if I load the blade in the edgewise direction, the difference in the torsional rotation at the blade tip becomes unphysical.

Could these differences be caused by issues 2 and 3 you described?

Regards, Jelmer

[andrew-platt]

The torsion deflection plot looks very similar to what I've seen with issues 2 and 3 described above (most likely a combination of the two). The ripple at the tip shows characteristics of the shape function for the last finite element node in the blade. That could either be caused either by incorrect twist associated with the mass and stiffness matrices near the tip from the cubic spline (issue 3 above), or from the incorrect reference curvature (issue 3 above -- though this is more noticeable with non-zero azimuth).

I haven't had a chance to get back to working on BeamDyn yet due to precedence of other projects. I'll reply when I have something further on this.

[JelmerPolman]

Hi Andy,

thank you for confirming.

Regards, Jelmer

[Leo-Antunes]

Dear all,

I am facing the same type of problem exposed by Jelmer, however using OpenFAST V2.3. Both codes from master and dev branch give the "solution does not converge" problem with BeamDyn. The time step 0.0002 works fine with the element order of 5, but the simulation takes too much time using this configuration. I would like to know with this is a normal time step or there is something I can do to improve the BeamDyn convergence.

Best regards,

Leonardo

[jjonkman]

Dear @Leo-Antunes,

Does the blade geometry you've defined via the key points in BeamDyn include built-in precurve or presweep? If so, is the reference curve smooth, e.g., fit with a low-order polynomial? If the reference curve is not smooth, then a high curvature-related stiffness could be introduced in your model, requiring a small time step.

Are all of the mass and stiffness entries set to physically meaningful values? If not, perhaps a high frequency is being introduced in your model, requiring a small time step.

Have you used the linearization capability of OpenFAST to identify the natural frequencies inherent in your model, including BeamDyn blade frequencies?

Best regards,

[Leo-Antunes]

Dear @jjonkman,

The blade doesn't have precurve or presweep and the mass and stiffness matrices seem to have reasonable values. About the linearization, I am having some difficulties to get the natural frequencies since it is a two-bladed model, I am working on it, but I am not familiar with Floquet theory yet,

However, I noted that reducing the damping coefficient the solution converges with higher values of time step. The initial values of the damping coefficients were 0.01s and time step for converging was 0.0003 s (DOF: Blades and yaw, 21 RPM) and element order of 7 when I changed the damping for 1E-03, it converges with DT = 0.00035s, using 1E-04s for damping coefficients the time step of 0.0009s is enough and with damping coefficients of 1E-05 s, even with DT = 0.02s, the simulation works fine.

Replacing my blade input file for the NREL 5MW 's file, the convergence was reached only with a DT of 0.0002s for the conditions quoted above. That's why I suppose it is not a problem with the mass and stiffness entries.

In the first post, I mentioned the time step of 0.0002s for convergence using an element order of 5, but these parameters were for a model with the tower's DOF set true.

Best Regards,

Leonardo

[jjonkman]

Dear @Leo-Antunes,

Regarding the linearization analysis, for the purposes of estimated the time-step requirements, you can likely get by with setting the rotor speed to zero, meaning the rotor is parked/idling and Floquet theory is not required to perform the eigenanalysis.

I'm not sure why the time step seems very sensitive to the stiffness-proportional damping coefficient.

Best regards,

[Leo-Antunes]

Dear @jjonkman,

When I try to perform the linearization the 0 RPM and stedy wind the following error appears:

Best regards,

[jjonkman]

Dear @Leo-Antunes,

When linearizing a parked/idling rotor to calculate natural frequencies, you can disable aerodynamics (CompInflow = CompAero = 0). If you want to include aerodynamic drag with CompInflow = CompAero = 1 to calculate the damping, you should disable the wake calculation in AeroDyn with WakeMod = 0 when the rotor is parked/idling (the wake model is only valid for an operational rotor).

Best regards,

[Leo-Antunes]

Dear @jjonkman,

Running the linearization, the highest full system natural frequency is 25.244 Hz. Taking the rule of thumb you mentioned at NREL forum for time step, dt <= 1/(10*f_max), good value for dt would be 0.004s, am I right?

Best regards,

[jjonkman]

Correct.

[Leo-Antunes]

Dear @jjonkman,

Continuing to investigate the reasons for the "Solution does not converge" problem, I noticed that increasing the hub mass and nacelle inercia the solution converges, but the values has to be too larger (almost 10E4 times larger). The dependency of these factors were also notable in the NREL 5 MW turbine: when I decrease the hub mass and nacelle inercia, setting the values of the model I am working on it (900 kg for hub mass and 3200 kg.m^2 for nacelle) , the same error "Solution does not converge after the maximum number of iterations" has shown. In this simulation I left the yaw free.

Best Regards,

[jjonkman]

Dear @Leo-Antunes,

These results make sense to me. Increasing the hub and nacelle mass/inertia will reduce the full-system natural frequencies, allowing for a larger time step. Conversely, reducing the hub and nacelle mass/inertia will increase the full-system natural frequencies, requiring a smaller time step.

Best regards,