View Issue Details
| ID | Project | Category | View Status | Date Submitted | Last Update |
|---|---|---|---|---|---|
| 0003366 | OpenFOAM | Contribution | public | 2019-10-09 14:19 | 2020-02-05 20:12 |
| Reporter | handrake0724 | Assigned To | henry | ||
| Priority | normal | Severity | minor | Reproducibility | have not tried |
| Status | closed | Resolution | suspended | ||
| Product Version | dev | ||||
| Summary | 0003366: fifth order Runge-Kutta implementation of rigydBodySolver | ||||
| Description | I wrote a Runge-Kutta method for integrating rigid body motion. In my field like offshore engineering, fourth or fifth order Runge-Kutta method are usually used. I have waited for the method to be introduced in the main tree and implemented one of the method myself. Please review if it is useful. | ||||
| Tags | No tags attached. | ||||
|
|
|
|
|
The Newmark scheme was implemented specifically for ship and offshore simulations on request from a sponsor, can you provide a case or cases which demonstrate the benefit of Runge-Kutta? |
|
|
I am not a expert for schemes especially Newmark and knows superfical things like: Newmark is 2nd order and popular Runge-Kutta is 4/5 order accuracy. Newmark is usually popular in FEM field. Newmark is lighter than Runge-Kutta. Newmark is direct solver of 2nd order ODE but Runge-Kutta is 1st order ODEs in state space I know your sponsor was from DNV GL. I guess they are using Newmark for solving motion. my colleagues who are naval archtectures are prefer to use RK4/5 for solving motion. Also offshore engineers I am acquainted with in Houston are also using RK4/5. There maybe reasons I don't know yet but considering my surroundings, choice of rigid body motion solver for ship hydrodynamics seems like somewhat preference. Thus, I think it might be beneficial to have one more method in that it can broaden choice of motion solvers to users, especially in my field. |
|
|
Can you provide a test case which shows the advantage of using Runge-Kutta 4/5 or any other scheme for the rigid body integration? If there is sensitivity to the integration method there are many many other choices we could investigate for this purpose like lower-order longer integration step Runge-Kutta, semi-implicit approaches or maybe even interface to the ODE integration package in OpenFOAM: Euler EulerSI RKCK45 RKDP45 RKF45 rodas23 rodas34 Rosenbrock12 Rosenbrock23 Rosenbrock34 seulex SIBS Trapezoid I have not done this so far because there was not a clear need for it but if you can provide cases which show this need we can investigate what schemes are most appropriate for a wider range of application. |
|
|
Thanks for considering this seriously. I am not sure if I could find a test case in wave induced motion problem but I will try to find in other fields. In astrophysics simulation, I was told numerical integrator is of great importance when it comes to solving a body or bodies revoluting around a body. Maybe such cases could be a candiate. The other day, I tried to find a usage in ODE package but it seemed prepared for lagrangian simulation and It looked like not an easy job to use in Eulerian simulation out of the box. So, it would be useful if the ODEpackage would have interface for Eulerian simulation. |
|
|
> In astrophysics simulation, I was told numerical integrator is of great importance when it comes to solving a body or bodies revoluting around a body. The symplectic integrator is often used for this purpose as it simultaneously conserves linear and angular momentum which RK methods don't. Note that when used in conjunction with the OpenFOAM fluid solvers it is not clear that higher than 2nd order is appropriate given that the fluid side is integrated with 2nd order at most and the interaction with the rigid body is lagged. If higher than 2nd order is required the time integration of the fluid side should also be changed and the interaction between the fluid and the rigid body included within the time integration sub-steps to ensure the coupled system is integrated at the higher order. |
|
|
> The other day, I tried to find a usage in ODE package but it seemed prepared for lagrangian simulation No, the ODE package was originally developed for chemistry integration, it has never been used for Lagrangian tracking. It does however have a general interface and if you need implicit integration requires the Jacobian but for RK and other explicit sub-step methods the Jacobian is not needed. If lots of different ODE intergrators are needed for rigidBody it would make sense for it to be interfaced to the current ODE system to avoid unnecessary code duplication. Currently we have no case or sponsors who need this functionality but if you would like to contribute this interface we can test it and include it if there is a general need. |
|
|
Yes you are right. it was for chemistry. my memory was not good since it was quite ago. I will take a look at how to add interface to ODE packages. |
|
|
henry, I wrote a little for ode support in rigidBodySolver. I took a similar structure of chemistrySolver such that rigidBodyMotion is inherited from ODESystem as well as from rigidBodyModel. also added is ode rigidBodySolver class. Now that rigidBodyMotion is inherited from ODESystem, it have to implement derivatives and jacobian member functions. I am not sure what should be in there. It would be helpful with your advice. please review the attached file to check if I am in the right track. patch.diff (3,417 bytes)
diff --git a/src/rigidBodyDynamics/Make/files b/src/rigidBodyDynamics/Make/files
index fd778ff9a..addab8e62 100644
--- a/src/rigidBodyDynamics/Make/files
+++ b/src/rigidBodyDynamics/Make/files
@@ -52,5 +52,6 @@ rigidBodySolvers/rigidBodySolver/newRigidBodySolver.C
rigidBodySolvers/symplectic/symplectic.C
rigidBodySolvers/Newmark/Newmark.C
rigidBodySolvers/CrankNicolson/CrankNicolson.C
+rigidBodySolvers/ode/ode.C
LIB = $(FOAM_LIBBIN)/librigidBodyDynamics
diff --git a/src/rigidBodyDynamics/Make/options b/src/rigidBodyDynamics/Make/options
index 79be6f3a7..bd4c3eb7b 100644
--- a/src/rigidBodyDynamics/Make/options
+++ b/src/rigidBodyDynamics/Make/options
@@ -1,3 +1,5 @@
-EXE_INC =
+EXE_INC = \
+ -I$(LIB_SRC)/ODE/lnInclude
-LIB_LIBS =
+LIB_LIBS = \
+ -lODE
diff --git a/src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotion.C b/src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotion.C
index db086b1c6..c2619a3cf 100644
--- a/src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotion.C
+++ b/src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotion.C
@@ -45,6 +45,7 @@ void Foam::RBD::rigidBodyMotion::initialize()
Foam::RBD::rigidBodyMotion::rigidBodyMotion()
:
rigidBodyModel(),
+ ODESystem(),
motionState_(*this),
motionState0_(*this),
aRelax_(1.0),
@@ -304,4 +305,25 @@ Foam::tmp<Foam::pointField> Foam::RBD::rigidBodyMotion::transformPoints
}
+void Foam::RBD::rigidBodyMotion::derivatives
+(
+ const scalar time,
+ const scalarField& y,
+ scalarField& dydx
+) const
+{
+
+}
+
+void Foam::RBD::rigidBodyMotion::jacobian
+(
+ const scalar time,
+ const scalarField& y,
+ scalarField& dfdx,
+ scalarSquareMatrix& dfdy
+) const
+{
+
+}
+
// ************************************************************************* //
diff --git a/src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotion.H b/src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotion.H
index 8b2d04efb..2070e99d8 100644
--- a/src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotion.H
+++ b/src/rigidBodyDynamics/rigidBodyMotion/rigidBodyMotion.H
@@ -50,6 +50,7 @@ SourceFiles
#include "rigidBodyModelState.H"
#include "pointField.H"
#include "Switch.H"
+#include "ODESystem.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -67,7 +68,8 @@ class rigidBodySolver;
class rigidBodyMotion
:
- public rigidBodyModel
+ public rigidBodyModel,
+ public ODESystem
{
friend class rigidBodySolver;
@@ -146,6 +148,31 @@ public:
// given body
spatialTransform X00(const label bodyId) const;
+ //- Return the number of equations in the system
+ virtual label nEqns() const
+ {
+ // nEqns = velocity 3 components and motion 3 components
+ return 6;
+ }
+
+ //- Calculate the derivatives in dydx
+ virtual void derivatives
+ (
+ const scalar time,
+ const scalarField& y,
+ scalarField& dydx
+ ) const;
+
+ //- Calculate the Jacobian of the system
+ // Need by the stiff-system solvers
+ virtual void jacobian
+ (
+ const scalar time,
+ const scalarField& y,
+ scalarField& dfdx,
+ scalarSquareMatrix& dfdy
+ ) const;
+
// Edit
|
|
|
derivatives member function would be implemented as follows: void Foam::RBD::rigidBodyMotion::derivatives ( const scalar time, const scalarField& y, scalarField& dydx ) const { dydx[0] = motionState0_.qDdot()[0]; dydx[1] = motionState0_.qDdot()[1]; dydx[2] = motionState0_.qDdot()[2]; dydx[3] = motionState0_.qDot()[0]; dydx[4] = motionState0_.qDot()[1]; dydx[5] = motionState0_.qDot()[2]; } |
|
|
Before going to all this effort we need to have a purpose for the implementation and more importantly the maintenance effort. Currently we have no cases or sponsors who need this; the ODE solvers provided with rigidBody have proved sufficient. Will your company fund this development and maintenance of the resulting code? Can they provide test and validation cases for which this development is required? |
|
|
as you said, rigidBodyDynamics is a general purpose class. so it can be used in other purpose other than CFD. confining to CFD, current prepared solvers up to 2nd order accuracy might be enough as you said. as a naval architecture and offshore engineer, I found another usage of this class in maneuvring simulations where governing equation is based on a perturbation model instead of doing CFD. manoeuvring simulations are usually used for simulator for training purpose or simulating collision avoidance near habor. in that case, higher order numerical integrators are required to track the position, i.e. 4/5 RK is a defacto standard in this field. I found that people around me are considering OpenFOAM toolkit for (big) data processing tool or general purpose toolkit, in addition to CFD toolkit. IMO, in that sense it would be more general and flexible and have more applications if rigidBodyDynamics class incorporates ODE package. If you think this feature needs a lot of work and efforts beyond a simple patch, I could come up with a JDP, maybe late next year. In that case, it would be helpful to plan JDP work scope if you let me know approximate quotation. |
|
|
If higher-order is required then you will need to handle the integration with the fluid fields higher-order as well, how did you do that? Did you compare the results with the symplectic integrator? Can you provide the test case which demonstrates the need for 4/5 RK? JDP? |
|
|
manoeuvring simulation is not solving flow field. all the hydrodynamic loads are modeled in perturbation model. Thus, It only needs rigid body dynamics and initial state variables. maybe I can write a test code implementing one of manoeuvring models using manoeuvring test data in public domain. JDP is joint development project. |
|
|
The integration via the general ODE solver system provides other options like sub-stepping, error control, semi-implicit etc. some or all of which may be useful for rigid body dynamics but currently then interface to CFD is not designed to support these but it could be done if there is a demonstrated need and funding as it will be quite complex to implement and maintain. If the rigid body dynamics library is used for other purposes and further development is required this can be undertaken given a specification for the need and funding. I had an inefficient prototype implementation of contact constraints but this is REALLY complex to do efficiently. The original sponsor of the rigid body dynamics library chose not to provide funding for the further development and maintenance work that they had promised after they had what they needed for their particular project. |
|
|
Ok, thank you for letting me know your situations. Do you think it is possible to make the rigidbodyDynamics class efficiently and have an idea? About the ODE part, I think I need to study the ODE and rigidBodyDynamics classes and make an application to manoeuvring models. So I would like to hold the discussion about it for now. I am now involved in joint industry project (JIP) developing efficient numerical wave tank simulation for offshore field during next year. Me and coordinator of the JIP have a plan in mind to have the developed code integrated with OpenFOAM. So end of next year I might make a proposal for funding those works. In that case, I would be able to add the scope improving rigid body dynamics and if possible, ODE. to be clear, is the fund for maintenance work made annually? if it is, it might be not possible. |
|
|
> Do you think it is possible to make the rigidbodyDynamics class efficiently and have an idea? The current rigidBodyDynamics implementation is efficient. > to be clear, is the fund for maintenance work made annually? if it is, it might be not possible. https://openfoam.org/news/funding-2019/ https://openfoam.org/maintenance/ |
|
|
Pending funding and suitable test cases to demonstrate the need for this addition. |
| Date Modified | Username | Field | Change |
|---|---|---|---|
| 2019-10-09 14:19 | handrake0724 | New Issue | |
| 2019-10-09 14:19 | handrake0724 | File Added: RK5.tar.gz | |
| 2019-10-09 15:49 | henry | Note Added: 0010824 | |
| 2019-10-10 03:41 | handrake0724 | Note Added: 0010825 | |
| 2019-10-10 05:32 | henry | Note Added: 0010826 | |
| 2019-10-10 07:53 | handrake0724 | Note Added: 0010827 | |
| 2019-10-10 08:08 | henry | Note Added: 0010828 | |
| 2019-10-10 08:13 | henry | Note Added: 0010829 | |
| 2019-10-10 08:49 | handrake0724 | Note Added: 0010830 | |
| 2019-10-11 02:25 | handrake0724 | File Added: patch.diff | |
| 2019-10-11 02:25 | handrake0724 | Note Added: 0010832 | |
| 2019-10-11 03:46 | handrake0724 | Note Added: 0010833 | |
| 2019-10-11 07:27 | henry | Note Added: 0010834 | |
| 2019-10-11 19:45 | handrake0724 | Note Added: 0010835 | |
| 2019-10-11 20:03 | henry | Note Added: 0010836 | |
| 2019-10-11 20:12 | handrake0724 | Note Added: 0010837 | |
| 2019-10-11 20:20 | henry | Note Added: 0010838 | |
| 2019-10-11 21:10 | handrake0724 | Note Added: 0010839 | |
| 2019-10-11 22:22 | henry | Note Added: 0010840 | |
| 2020-02-05 20:12 | henry | Assigned To | => henry |
| 2020-02-05 20:12 | henry | Status | new => closed |
| 2020-02-05 20:12 | henry | Resolution | open => suspended |
| 2020-02-05 20:12 | henry | Note Added: 0011119 |
