#### View Issue Details

ID Project Category View Status Date Submitted Last Update 0002211 OpenFOAM Bug public 2016-08-23 17:47 2019-06-11 13:33 yuanchuanliu henry normal minor have not tried resolved fixed GNU/Linux Ubuntu 16.04 0002211: Negative weight before averaging in rigidBodyMeshMotion In the transformPoints function of the rigidBodyMotion class for multiple bodies, before the average function is called to produce an averaged septernion, the weights of a point regarding multiple bodies are "re-normalised" to produce an array of new weights. However, some of the new weights could be negative. I am not sure if a negative weight would still be valid for averaging. Let us do some math using the logic I understand from the code. For a system of three rigid bodies, the original weight List for any point then has a total of 4 elements (body number + 1), of which the first three weights for the point are assumed to be: 0.9, 0.8 and 0.3. The new weight list can be calculated as follows: 1. Sum 1-w: 0.1 + 0.2 + 0.7 = 1.0; 2. Find the maximum weight: 0.9; 3. Calculate lamda: (4 -1 - 0.9)/1.0 = 2.1; 4. "Re-normalise" weights:     a) 1 - 2.1*(1 - 0.9) = 0.79;     b) 1 - 2.1*(1 - 0.8) = 0.58;     c) 1 - 2.1*(1 - 0.3) = -0.47; (NEGATIVE) 5. Sum new weights: 0.79 + 0.58 -0.47 = 0.9 (equal to original maxw); 6. Calculate the last weight for the far-field: 1- 0.9 = 0.1. When one of the original weights is 1.0, whatever the other two weights originally are, the new weights would always be: 1.0 (this is well preserved), w and -w. If the corresponding three septernions are s1, s2 and s3, after averaging, would the point still moves along with the body producing a weight of 1.0 for it? I am raising this issue as I am modelling a three-body system connected by rotational joints and the mesh always look badly distorted after a few time steps (even the initial body shape cannot be preseved). I am not confident my settings in my constant/dynamicMeshDict are correct but this logic is something I also want to understand. No tags attached.

#### Activities

 2016-08-23 18:05 manager   ~0006727 What would you do to avoid negative weights? If you apply some "fix" to the negative weights is the resulting motion closer to expectation? 2016-08-23 18:15 reporter dynamicMeshDict.txt (2,630 bytes)    ```/*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | | \\ / O peration | Version: 4.x | | \\ / A nd | Web: www.OpenFOAM.org | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ FoamFile { version 2.0; format ascii; class dictionary; object dynamicMeshDict; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // dynamicFvMesh dynamicMotionSolverFvMesh; motionSolverLibs ("librigidBodyMeshMotion.so"); solver rigidBodyMotion; rigidBodyMotionCoeffs { report on; solver { type Newmark; } test true; nIter 2; accelerationRelaxation 0.4; bodies { body1 { type rigidBody; parent root; centreOfMass (0.5 0 0); mass 412.73; inertia (40 0 0 921 0 921); transform (1 0 0 0 1 0 0 0 1) (0.5 0 0); joint { type composite; joints ( { type Py; } { type Rz; } ); } patches (body1); innerDistance 0.1; outerDistance 1; } body2 { \$body1; parent body1; centreOfMass (1.5 0 0); transform (1 0 0 0 1 0 0 0 1) (1 0 0); patches (body2); } body3 { \$body1; parent body2; centreOfMass (2.5 0 0); transform (1 0 0 0 1 0 0 0 1) (1 0 0); patches (body3); } } restraints { spring1 { type linearSpring; body body1; anchor (0 0.5 0); refAttachmentPt (0 0 0); stiffness 1e5; damping 0; restLength 0.4; } } } // ************************************************************************* // ``` dynamicMeshDict.txt (2,630 bytes) 2016-08-23 18:16 reporter   ~0006728 Last edited: 2016-08-23 18:23 Do you think negative weights would be OK? To be honest, I have not fully understood how to properly set up the dynamicMeshDict for a three-body system. If you can help me with the dictionary I uploaded, I can use it to try and see if it can be improved. Cheers. Edit: my model consists of three bodies connected by rotational joints. The spring applied on the first body is used to make it move as well as two other bodies. 2016-08-23 18:31 manager   ~0006729 The current method does not preclude negative weights and it is not clear if it should be changed to avoid them and if so by what means. If you can provide an alternative method which avoids negative weights and improves the motion I would be happy to include it in OpenFOAM-dev. 2016-08-24 15:51 reporter   ~0006744 I made several attempts but none of them could entirely improve the situation. The best I can get till now is to maintain the quality of the mesh surrouding the bodies, but the biggest problem of not being able to preseve the initial body shape is still around. Please see the attached images showing the mesh around the bodies. zero.png shows the initial body shape and mesh. two-orig.png shows the mesh at T = 2 s using the algorithm provided with the release. two.png shows the mesh using the currently tested method. The present method first applys the slerp function on the septernion related to each body with the corresponding weight to get a new septernion. For three bodies, there are three new septernions for a given point. Then apply the average function on the new septernions while the new weights are obtained by normalising the original weights with their sum. I cannot think of a way to preserve the initial body shape. Forcing the points within the innerDistance of a body to move with it will only make it worse. 2016-08-24 15:52 reporter zero.png (105,939 bytes)    zero.png (105,939 bytes) 2016-08-24 15:52 reporter two-orig.png (140,370 bytes)    two-orig.png (140,370 bytes) 2016-08-24 15:52 reporter two.png (144,152 bytes)    two.png (144,152 bytes) 2016-08-24 16:00 manager   ~0006747 > Forcing the points within the innerDistance of a body to move with it will only make it worse. Why? The purpose of the innerDistance and outerDistance is to ensure that the motion of neighbouring bodies to not cause distortion of the initial body shapes. This approach has worked for the cases I have tested. 2016-08-24 16:38 reporter   ~0006749 I totally understand the purpose of defining the innerDistance and outerDistance(points within the innerDistance will move with the body while those out of outerDistance will stay stationary). It is just that the initial body shape is never preseved for my multi-body cases (the slerp approach works for single body cases of course) even when I used the algorithm provided in the official release. Maybe my dynamicMeshDict settings for the multibody system are not correct? It would be very helpful if you can provide a simple case you successfully tested. If the current method works well for multibody problems, negative weights would not be a issue at all. 2016-08-24 16:59 manager rigidBodySHM2.tgz (3,123 bytes) 2016-08-24 16:59 manager   ~0006750 I have attached a simple 2-body 3D test case. 2016-08-24 17:46 reporter   ~0006751 Thank you for providing a test case. I just realised that a two-body system would be OK but it is not true for a system with more than two bodies. Please see my 3-body case attached. If we do the same math again. For a 2-body system, the first two weights for the point are assumed to be: 1 and 0.8. The new weight list can be calculated as follows: 1. Sum 1-w: 0 + 0.2 = 0.2; 2. Find the maximum weight: 1; 3. Calculate lamda: (3 - 1 - 1)/0.2 = 5; 4. "Re-normalise" weights:     a) 1 - 5*(1 - 1) = 1;     b) 1 - 5*(1 - 0.8) = 0; 5. Sum new weights: 1 + 0 = 1 (equal to original maxw); 6. Calculate the last weight for the far-field: 1 - 1 = 0. In fact, for a 2-body system, when the original weight is 1 for a point, whatever the other weight is originally, it will always be 0 after "re-normalisation", meaning the points within innerDistance will indeed move with the body as expected. This is due to the rule set beforehand that sum of new weights will always be equal to the max of the original weight list. For a system with more than two bodies, however, it will not always hold true as has been stated in Steps To Reproduce. 2016-08-24 17:46 reporter RBD-test.zip (603,442 bytes) 2016-08-24 17:58 manager   ~0006752 I used the following reference for the averaging of quaternions: http://malcolmdshuster.com/FC_MarkleyMortari_Girdwood_1999_AAS.pdf but I only implemented the simplest method; the more complex methods might be better for many bodies. 2016-08-24 18:09 reporter   ~0006753 Thank you very much for providing a reference paper. It is definitely worth reading it carefully. Please close the issue for now if you wish. Once I have something new to share with you, I will reopen it. Cheers. 2016-09-14 10:02 manager RBD-laplacian.tgz (603,683 bytes) 2016-09-14 10:06 manager   ~0006860 New in OpenFOAM-dev: commit 2c90bd2ee6d674c4f693d8e9294cceca8f311d98 Author: Henry Weller Date: Wed Sep 14 09:59:02 2016 +0100     rigidBodyMeshMotionSolver: experimental nDoF mesh-motion solver supporting the displacement-based elliptic solvers          Specification for the tutorials/multiphase/interDyMFoam/ras/floatingObject case:          dynamicFvMesh dynamicMotionSolverFvMesh;          motionSolverLibs ("librigidBodyMeshMotion.so" "libfvMotionSolvers.so");          solver rigidBodyMotionSolver;          rigidBodyMotionSolverCoeffs     {         report on;              meshSolver         {             solver displacementLaplacian;                  displacementLaplacianCoeffs             {                 diffusivity inverseDistance (floatingObject);             }         }     .     .     . I have also attached a version of your test-case using the laplacian motion solver: RBD-laplacian.tgz However, it would still be good to resolve the issues with the SLERP interpolation for >2 bodies. 2016-09-14 14:16 manager   ~0006861 I have implemented a more general transform averaging procedure: commit 47811eae8e9b661250cbc0ea151a89debcc9783e could you test it on your cases and report back? 2016-09-14 17:08 reporter   ~0006862 After some test, I can confirm that the updated implementation is able to preserve the rigid body shape and also points inside the innerDistance will move with the corresponding body as expected although the values for the innerDistance and outerDistance might need some trials. Thanks for the update. 2016-09-14 17:11 manager   ~0006863 Have you also tried with the laplacian motion solver? 2016-09-15 14:20 reporter   ~0006864 Sorry for my late response. Yes, I also tried the new laplacian mesh motion solver and it worked well. And I love your brilliant idea of including the meshSolver as a pointer in the rigidBodyMeshMotionSolver class so that different solvers might be selected to solve the mesh motion. Cheers. 2016-09-15 14:37 manager   ~0006865 Resolved in OpenFOAM-dev by commit 47811eae8e9b661250cbc0ea151a89debcc9783e Resolved in OpenFOAM-4.x by commit de29a3309ecea18ac159d710ff3fa346850d15d1

#### Issue History

2016-08-23 17:47 yuanchuanliu New Issue
2016-08-23 18:05 henry Note Added: 0006727
2016-08-23 18:15 yuanchuanliu File Added: dynamicMeshDict.txt
2016-08-23 18:16 yuanchuanliu Note Added: 0006728
2016-08-23 18:23 yuanchuanliu Note Edited: 0006728
2016-08-23 18:31 henry Note Added: 0006729
2016-08-24 15:51 yuanchuanliu Note Added: 0006744
2016-08-24 15:52 yuanchuanliu File Added: zero.png
2016-08-24 15:52 yuanchuanliu File Added: two-orig.png
2016-08-24 15:52 yuanchuanliu File Added: two.png
2016-08-24 16:00 henry Note Added: 0006747
2016-08-24 16:38 yuanchuanliu Note Added: 0006749
2016-08-24 16:59 henry File Added: rigidBodySHM2.tgz
2016-08-24 16:59 henry Note Added: 0006750
2016-08-24 17:46 yuanchuanliu Note Added: 0006751
2016-08-24 17:46 yuanchuanliu File Added: RBD-test.zip
2016-08-24 17:58 henry Note Added: 0006752
2016-08-24 18:09 yuanchuanliu Note Added: 0006753
2016-09-14 10:02 henry File Added: RBD-laplacian.tgz
2016-09-14 10:06 henry Note Added: 0006860
2016-09-14 14:16 henry Note Added: 0006861
2016-09-14 17:08 yuanchuanliu Note Added: 0006862
2016-09-14 17:11 henry Note Added: 0006863
2016-09-15 14:20 yuanchuanliu Note Added: 0006864
2016-09-15 14:37 henry Assigned To => henry
2016-09-15 14:37 henry Status new => resolved
2016-09-15 14:37 henry Resolution open => fixed
2016-09-15 14:37 henry Fixed in Version => 4.x
2016-09-15 14:37 henry Note Added: 0006865