View Issue Details
ID  Project  Category  View Status  Date Submitted  Last Update 

0002211  OpenFOAM  Bug  public  20160823 17:47  20190611 13:33 
Reporter  yuanchuanliu  Assigned To  henry  
Priority  normal  Severity  minor  Reproducibility  have not tried 
Status  resolved  Resolution  fixed  
Platform  GNU/Linux  OS  Ubuntu  OS Version  16.04 
Summary  0002211: Negative weight before averaging in rigidBodyMeshMotion  
Description  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 "renormalised" 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.  
Steps To Reproduce  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 1w: 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. "Renormalise" 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 farfield: 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?  
Additional Information  I am raising this issue as I am modelling a threebody 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.  
Tags  No tags attached.  

What would you do to avoid negative weights? If you apply some "fix" to the negative weights is the resulting motion closer to expectation? 

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; } } } // ************************************************************************* // 

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 threebody 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. 

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 OpenFOAMdev. 

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. twoorig.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. 







> 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. 

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 multibody 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. 



I have attached a simple 2body 3D test case. 

Thank you for providing a test case. I just realised that a twobody system would be OK but it is not true for a system with more than two bodies. Please see my 3body case attached. If we do the same math again. For a 2body 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 1w: 0 + 0.2 = 0.2; 2. Find the maximum weight: 1; 3. Calculate lamda: (3  1  1)/0.2 = 5; 4. "Renormalise" 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 farfield: 1  1 = 0. In fact, for a 2body system, when the original weight is 1 for a point, whatever the other weight is originally, it will always be 0 after "renormalisation", 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. 



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. 

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. 



New in OpenFOAMdev: commit 2c90bd2ee6d674c4f693d8e9294cceca8f311d98 Author: Henry Weller <http://cfd.direct> Date: Wed Sep 14 09:59:02 2016 +0100 rigidBodyMeshMotionSolver: experimental nDoF meshmotion solver supporting the displacementbased 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 testcase using the laplacian motion solver: RBDlaplacian.tgz However, it would still be good to resolve the issues with the SLERP interpolation for >2 bodies. 

I have implemented a more general transform averaging procedure: commit 47811eae8e9b661250cbc0ea151a89debcc9783e could you test it on your cases and report back? 

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. 

Have you also tried with the laplacian motion solver? 

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. 

Resolved in OpenFOAMdev by commit 47811eae8e9b661250cbc0ea151a89debcc9783e Resolved in OpenFOAM4.x by commit de29a3309ecea18ac159d710ff3fa346850d15d1 
Date Modified  Username  Field  Change 

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