OpenFAST / openfast

Main repository for the NREL-supported OpenFAST whole-turbine and FAST.Farm wind farm simulation codes.
http://openfast.readthedocs.io
Apache License 2.0
673 stars 455 forks source link

Potential for Under Integration Using Trapezoidal Quadrature in BeamDyn #183

Open jjonkman opened 5 years ago

jjonkman commented 5 years ago

BeamDyn does not currently check for consistency between the number of trapezoidal quadrature points and the number of finite-element nodes. If there are too few trapezoidal quadrature points for the number of finite-element nodes, the spatial integrals will be under integrated and BeamDyn will produce incorrect solutions or will fatally crash with a singular matrix error. A check should be added within BeamDyn to ensure that the use of trapezoidal quadrature does not lead to under integration.

There is no problem with Gaussian quadrature.

andrew-platt commented 5 years ago

I implemented a simple set of test criteria using the shape functions (stored in p%Shp) in the Envision version of BeamDyn that appears to catch most of these issues. There is undoubtedly a more robust mathematical method for testing this.

The following criteria were found through experimentation with three blade designs and approximately 90 test cases with various combinations of order of element and refinement.

  1. if the sum of QP contributions to an FE is negative. This is guaranteed to fail. sum(p%Shp(FE,:)) < 0.0
  2. if the normalized sum of the QP contributions to an FE is small (normalize by total FE nodes). This may run but is likely to diverge sum(p%Shp(FE,:)) < 0.020
  3. if the ratio of min and max FE QP contributions is very low. This may run and return reasonable approximations, but will diverge under large loads (>10 g or high aero loading conditions). This value was found also through observation and has at present no theoretical backing. I suspect the reason this test works is that it indicates that a given FE node is proportionally under-represented. The limit is found empirically. sum(p%Shp(FE,:)) / maxval(sum(p%Shp(:,:),2)) < 0.13