View Issue Details

Jump to NotesJump to History
IDProjectCategoryView StatusDate SubmittedLast Update
0002211OpenFOAMBugpublic2016-08-23 17:472019-06-11 13:33
Reporteryuanchuanliu Assigned Tohenry  
PrioritynormalSeverityminorReproducibilityhave not tried
Status resolvedResolutionfixed 
PlatformGNU/LinuxOSUbuntuOS Version16.04
Summary0002211: Negative weight before averaging in rigidBodyMeshMotion
DescriptionIn 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.
Steps To ReproduceLet 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?
Additional InformationI 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.
TagsNo tags attached.

Activities

henry

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?

yuanchuanliu

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)   

yuanchuanliu

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.

henry

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.

yuanchuanliu

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.

yuanchuanliu

2016-08-24 15:52

reporter  

zero.png (105,939 bytes)   
zero.png (105,939 bytes)   

yuanchuanliu

2016-08-24 15:52

reporter  

two-orig.png (140,370 bytes)   
two-orig.png (140,370 bytes)   

yuanchuanliu

2016-08-24 15:52

reporter  

two.png (144,152 bytes)   
two.png (144,152 bytes)   

henry

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.

yuanchuanliu

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.

henry

2016-08-24 16:59

manager  

rigidBodySHM2.tgz (3,123 bytes)

henry

2016-08-24 16:59

manager   ~0006750

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

yuanchuanliu

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.

yuanchuanliu

2016-08-24 17:46

reporter  

RBD-test.zip (603,442 bytes)

henry

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.

yuanchuanliu

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.

henry

2016-09-14 10:02

manager  

RBD-laplacian.tgz (603,683 bytes)

henry

2016-09-14 10:06

manager   ~0006860

New in OpenFOAM-dev:

commit 2c90bd2ee6d674c4f693d8e9294cceca8f311d98
Author: Henry Weller <http://cfd.direct>
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.

henry

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?

yuanchuanliu

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.

henry

2016-09-14 17:11

manager   ~0006863

Have you also tried with the laplacian motion solver?

yuanchuanliu

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.

henry

2016-09-15 14:37

manager   ~0006865

Resolved in OpenFOAM-dev by commit 47811eae8e9b661250cbc0ea151a89debcc9783e
Resolved in OpenFOAM-4.x by commit de29a3309ecea18ac159d710ff3fa346850d15d1

Issue History

Date Modified Username Field Change
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