Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

SD: Implementing directional cosine matrices and section properties for rectangular members #1413

Merged
merged 21 commits into from Apr 14, 2023

Conversation

samuel-ramsahoye
Copy link
Contributor

@samuel-ramsahoye samuel-ramsahoye commented Jan 18, 2023

This PR is ready to be merged.

Compilation has only been experimented with via compilation of the SubDyn driver ("make subdyn_driver"). The physical correctness will be checked soon via some canonical tests.

No new tests have been added yet for this new functionality - existing tests may fail.

Feature or improvement description
SubDyn primary input files allow the user to specify "X" section properties and member cosines as tables but currently these inputs are not used. This PR aims to correct that.

This update should be backwards compatible with existing SubDyn primary input files. If the COSMID column is missing from the member specification, SubDyn's default member cosine convention is used. Moreover, the user can partially specify member cosines by specifying -1 in place of COSMID for the relevant rows: -1 is the fill value for the entire column if it's left unspecified.

"X" section properties for a member are used if MType=4. Whilst MType=4 members can be refined, the property ID at either joint must be the same. If not, a fatal exception occurs and an appropriate message is displayed. The definition of "X" section properties hasn't changed from previous implementations.

User-specified cosines are a direct replacement for the GetDirCos call in SD_FEM.f90 and so they are in fact the transpose of what would normally be considered the member cosine matrix. At the moment, no checks are in place to ensure the correctness of user-defined cosines.

Related issue, if one exists
#1404

Impacted areas of the software
This PR relates to the SubDyn module only.

Changes predominantly relate to SubDyn.f90, SD_FEM.f90 and FEM.f90.

Additional supporting information
The "./subdyn_cosine_arbitrary_example" has been used during initial testing. This code has only be experimented with via the SubDyn driver and not via the OpenFAST glue-code and compilation has only been checked via "make subdyn_driver".

Test results, if applicable
Tests will be added to this draft PR is due course.

@samuel-ramsahoye samuel-ramsahoye marked this pull request as draft January 18, 2023 14:21
@samuel-ramsahoye
Copy link
Contributor Author

@ebranlard Please see the draft PR for the member cosines and X section properties issue.

@samuel-ramsahoye
Copy link
Contributor Author

samuel-ramsahoye commented Jan 20, 2023

@bjonkman Thanks very much Bonnie!
@ebranlard I’ll put together some tests for this new functionality within a week.

@samuel-ramsahoye
Copy link
Contributor Author

@ebranlard I've found a couple of issues with my PR in my regression tests, the main one being that I don't think the stiffness matrices are actually being evaluated for members with X sections yet. I know what I need to modify in the code to fix this.

For the testing, to start with what I'm doing is taking existing SubDyn models and re-writing them by specifying member cosines and properties via X sections and checking the results are still the same.

@samuel-ramsahoye
Copy link
Contributor Author

samuel-ramsahoye commented Jan 20, 2023

@ebranlard There's an issue with the user-specified cosines in this PR. I set up a simple angled cantilever example, the modes of which should be the same irrespective of whether or not the member cosines are specified by the user or calculated by GetDirCos. I populated the COSM information by writing it out from SubDyn. The modal frequencies are the same but the displacements are different, except for the first mode. The relevant input files for this case are in the ./f_subdyn_test/02 folder if you want to take a look. You can reconstruct this behaviour by toggling the COSMID value of the single member between -1 and 1.

Copy link
Contributor

@ebranlard ebranlard left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for doing that, the changes looks really good. I mostly have a couple of comments regarding some conventions that I would personally prefer.

One thing I couldn't comment on but might be important:

  • I would recommend setting COSMs (in the registry) as R8Ki and not ReKi. Maybe this is the reason why you see different results
    Also:
  • I would recommend using another test case, where you keep the regular Beam type (3) and specify the directional cosine. (to avoid testing both the directional cosine feature and the anisotropic member feature).

Some important points to think about (and probably discuss with @jjonkman):

  • Do we agree that "anisotropic" (Aniso for short) is better than "X" ?
  • do we want to allow for division of anisotropic members?
  • Should we average the properties on both ends of the anisotropic member? (We need to think again of what is done in SubDyn when there are different properties on both node ends).
  • We need to make sure the convention for the directional cosine inputs is well documented (given the internal transpose required).

@@ -0,0 +1,26 @@
SubDyn Driver file for stand-alone applications
Copy link
Contributor

@ebranlard ebranlard Jan 20, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

NOTE: we will have to remove these input files as they should be placed in the r-test repository. It's fine to keep it here for now. I will squash this pull request at the end so the adding and deleting of the files will not appear in the list of commits.

My own preferences: rename the driver with extension .dvr.
I would personnaly not put the .yaml and .json files in the r-test when the time comes (I know it's present in other tests..).

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes I completely agree - currently just there incidentally during development.

modules/subdyn/src/SD_FEM.f90 Outdated Show resolved Hide resolved
@@ -58,6 +59,7 @@ MODULE SD_FEM
INTEGER(IntKi), PARAMETER :: idMemberBeam = 1
INTEGER(IntKi), PARAMETER :: idMemberCable = 2
INTEGER(IntKi), PARAMETER :: idMemberRigid = 3
INTEGER(IntKi), PARAMETER :: idMemberX = 4
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I know the input file use NXPropSets but I find this "X" slightly confusing. Maybe we replace "X" with "Aniso" (for anisotropic), or anything that would make sense. Then we'll likely change the input file to have NAnisoPropSets. You can hold off for now, I'll need @jjonkman's opinion.

modules/subdyn/src/SD_FEM.f90 Outdated Show resolved Hide resolved
modules/subdyn/src/SD_FEM.f90 Outdated Show resolved Hide resolved
modules/subdyn/src/SubDyn.f90 Outdated Show resolved Hide resolved
modules/subdyn/src/SubDyn.f90 Outdated Show resolved Hide resolved
modules/subdyn/src/SubDyn.f90 Outdated Show resolved Hide resolved
modules/subdyn/src/SubDyn.f90 Outdated Show resolved Hide resolved
@@ -0,0 +1,2 @@
1. vertical cantilever
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In the r-test, use README.md for these kind of files if you can.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes of course, will formalise with README files once I’m completely satisfied with the results

@samuel-ramsahoye
Copy link
Contributor Author

@ebranlard Thanks very much for your comments. I’ve made the majority of the changes you suggested. I believe the outstanding ones are conditional on @jjonkman’s opinion.

Your suggestion of increasing the number of bits used to store the direction cosine information worked when additionally used in conjunction with writing and reading the member cosines to a much higher precision. The precision with which the user needs to specify these cosines is very high in order for the results to match “perfectly”. This makes me think that specification via a rotational angle about the direction vector, rather than 9 high precision elements of a matrix, might be better in retrospect.

To summarise, results match when both:

  1. The user directly specifies CHS properties via equivalent X section properties
  2. The user specifies member direction cosines rather than using the internal evaluation in SubDyn

The next few steps for me (aside from any further recommendations from yourself):

  1. Formalise a set of canonical tests that separate the checking of each of these two new features independently as well as together
  2. I’d like to review the changes I’ve made again now all the functionality appears to be working as expected

Cheers,
Sam

@ebranlard
Copy link
Contributor

Thanks a lot for doing the changes! You'll see that the tests are failing in single precision, because FINLOCI (in FEM.f90) was written for ReKi or Int. The easiest is probably to add another interface to FINDLOCI (for SiKi) and change the existing one from ReKi to R8Ki.

All rotations in OpenFAST are stored as double precision as these parameters can indeed be quite sensitive. We could input the directional cosines as 3 Euler angles or 4 Euler parameters in the input file (preferable to only 3 parameters), but I would expect those to also have to be given with high precision.

@samuel-ramsahoye
Copy link
Contributor Author

@ebranlard For our particular use case, it's fine to assume the direction cosine of the local element Z-axis is the element direction vector. The orientation of the local element X and Y axes can then be determined by a single rotational angle about the local element Z-axis. I think this might have been what you originally suggested. This could be a member (and element) property rather than (a column in the input file in the member table) and we could then get rid of the direction cosines table.

We could probably do with some error handling on user-specified cosines: I'll try to add in the next few days.

Cheers,
Sam

@ebranlard
Copy link
Contributor

Hi @samuel-ramsahoye , yes, the rotation around the z axis was indeed what I had in mind for ease of input by the user. The user still have to think about/figure out where the "x" axis is placed by default by SubDyn. I would probably hold on the implementation for now until we take a decision, or at least not remove the direction cosine table.

@samuel-ramsahoye
Copy link
Contributor Author

@ebranlard Yes I’ll hold off. Does NREL provide a means to visualise OpenFAST models? It would be nice if the user could visually inspect the orientation of the local-element axes in SubDyn, amongst other things from other modules too. If not I can write something in Python.

@ebranlard
Copy link
Contributor

I hope to get back to you in a week or so for the remaining decisions.

Regarding visualization of the orientation, I've never done it, but my colleague @andrew-platt has. The only available option at the moment is to output the meshes to VTK using WrVTK and VTK_Type = 2 (for instance). Then in ParaView, you can use 3D glyph arrows to visualize each component of the orientation matrix (which is basically your frame vectors):
Screen Shot 2023-01-25 at 10 46 26 AM (1)

@andrew-platt
Copy link
Collaborator

Following on what @ebranlard mentioned, here is an example of what visualizing the orientations can look like (X = red, Y = yellow, Z = green). This example is from visualization of the blade root orientations from two simulations (lighter arrows are the reference simulation).

Screen Shot 2023-01-25 at 10 56 35 AM

@ebranlard
Copy link
Contributor

Hi @Sam-Ramsahoye

I think we assumed that the elements would mostly be rectangular, but that indeed that doesn't have to be the case. Then NonCirc" or "Arb" (for arbitrary) would be better.

Sorry for the confusion about isotropic. I was reading the input file of SubDyn and seeing that the "X-member" section was intended for "isotropic" material, which I believe is isotropic in the sense that the material is the same throughout the member (same density, E, and G), but indeed the geometrical properties are not axisymmetric and can vary along the span of the member.

Let me know if that makes more sense!
Thanks again for all your work,

Emmanuel

@Sam-Ramsahoye
Copy link

@ebranlard No worries, thanks for clarifying.

I can look at making some of those changes over the next couple of days. The only thing I might not have time to do is write in the interpolation of NonCirc properties. Currently, those elements are subdivided because the simplification that the start and end properties must be the same makes it simple.

@ebranlard
Copy link
Contributor

Great, thanks a lot for that, it's fine if you don't have the time to do the interpolation.

@verlivkra
Copy link

verlivkra commented Mar 8, 2023

@ebranlard Yes I’ll hold off. Does NREL provide a means to visualise OpenFAST models? It would be nice if the user could visually inspect the orientation of the local-element axes in SubDyn, amongst other things from other modules too. If not I can write something in Python.

@samuel-ramsahoye Did you find a way to visualize the member local coordinate system? I thought it was not possible to visualize SubDyn-members in Paraview as of now, but it would be very useful - both to visualize the SubDyn-members and see their local coordinate systems. Otherwise, maybe this could be useful to include in viz3Danim, @ebranlard?

@ebranlard
Copy link
Contributor

Hi @verlivkra, did the method mentioned by @andrew-platt not work for you? (I have not tried it..)

@ebranlard
Copy link
Contributor

Hi @Sam-Ramsahoye, have you made any progress on the remaining implementation detail? We would like to merge this pull request in the next couple of weeks. Thanks!

@Sam-Ramsahoye
Copy link

@ebranlard Apologies for my delay in getting back to you. I've been very busy recently but planning to look at this early next week.

Apologies for the delay again.

Sam

@verlivkra
Copy link

@ebranlard as far as I can see, no SubDyn-members are visualized when generating vtk-files from OpenFAST. I don't see how 3D-glyphs would then be available for the SubDyn-members. What I saw from @andrew-platt 's comment was only an example on how this could look for the blades, which I assume are not SubDyn-members. Did I miss something here?
If there is a better place to continue this discussion, please let me know.

@jjonkman
Copy link
Collaborator

Dear @verlivkra,

FYI: The issue preventing visualization of the substructure in ParaView (#776) has not been fixed yet. That said, you can visualize the Guyan and Craig-Bampton modes of the SubDyn model through the following website created by Emmanuel Branlard: https://ebranlard.github.io/viz3Danim/; to use this feature, simply load the .json summary file generated by SubDyn.

Best regards,

@ebranlard
Copy link
Contributor

Hi @samuel-ramsahoye
Do you have time to wrap up this pull request? We are hoping to put that in our 3.5.0 release in the coming weeks. If you don't have time, can you grant me access to your branch? (add me as a collaborator)
Could you also add a zip file in a post here with the small test cases you had put together?
Thanks!

@samuel-ramsahoye
Copy link
Contributor Author

Hi @ebranlard,

Really sorry I’ve been so slow. Unfortunately, I haven’t had time to have a look at this recently and a jacket project I had planned on studying in OpenFAST, where the functionality would have been useful, got postponed until May.

I’ll add you as a contributor and try and gather the (very) limited test cases I’ve put together.

Cheers,
Sam

@andrew-platt andrew-platt added this to the v3.5.0 milestone Apr 10, 2023
@ebranlard ebranlard marked this pull request as ready for review April 14, 2023 21:47
@ebranlard ebranlard self-requested a review April 14, 2023 22:48
modules/subdyn/src/SD_FEM.f90 Outdated Show resolved Hide resolved
modules/subdyn/src/SD_FEM.f90 Show resolved Hide resolved
modules/subdyn/src/SD_FEM.f90 Outdated Show resolved Hide resolved
modules/subdyn/src/SD_FEM.f90 Outdated Show resolved Hide resolved
@ebranlard ebranlard changed the title Implementing directional cosine matrices and section properties for rectangular members SD: Implementing directional cosine matrices and section properties for rectangular members Apr 14, 2023
@ebranlard ebranlard merged commit e5df74b into OpenFAST:dev Apr 14, 2023
29 checks passed
@Sam-Ramsahoye
Copy link

@ebranlard Apologies I didn't get around to putting together the tests for you, I've been very busy.

I've got a couple of uses for this functionality coming up so I'm going to collate the tests then - some additional checks that I've been performing seem to compare well.

Would it be cleaner to compare expected and calculated stiffness matrices rather than time-series outputs from quasi-static sims?

I'm aware this has now been merged but I've found a small bug in this update. If you don't specify any beams with regular circular sections, an error is thrown. This is because currently when the p%ElemProps is set for eType==idMemberBeamArb, D1 and D2 are referenced but undefined. In any case, these are erroneous references and should probably replaced with A**0.5 or some other equivalent diameter.

@ebranlard
Copy link
Contributor

Thanks for looking into it! Hopefully #1531 fixes this issue.

@andrew-platt andrew-platt mentioned this pull request May 23, 2023
19 tasks
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

7 participants