[Ship Designer]

By definition and by way of how the forces function object calculates forces and moments, these two vectors should be perpendicular to each other, i.e. their scalar product should be zero. I have found that in many of the cases I work on and some tutorials as well, that the total force and moment vectors output by the function object are not perpendicular. Their scalar product is not zero, with angles ranging from 53° to 91° for the cases I checked. This is bothering me because the discrepancy seems to be substantial and I need reasonably accurate moments for my work. The DTC tutorial case (fixed hull resistance) for example has an angle of 53° between F and M, see below as calculated in OF v8 at last time step:

forces forces write:
    sum of forces:
        pressure : (-2.66933 1748.01 3914.97)
        viscous  : (-12.2637 -0.197327 -0.279397)
        porous   : (0 0 0)
    sum of moments:
        pressure : (-532.767 29.4779 444.697)
        viscous  : (0.0187431 2.99189 -3.72562)
        porous   : (0 0 0)

The pressure force-moment scalar product is 1793925 with an angle of 53° and the viscous one is 0.221 with an angle of 89.8°, which is acceptable for the latter I think.

The forces function object calculates the moment by calculating the cross product d ⨯ F, where d is an array of the face centres corrected by the CofR dictionary entry, if present, and F the forces array. If I'm not mistaken this happens on line 382 of the forces.C file of OF v8:

void Foam::functionObjects::forces::applyBins
(
    const vectorField& Md,
    const vectorField& fN,
    const vectorField& fT,
    const vectorField& fP,
    const vectorField& d
)
{
    if (nBin_ == 1)
    {
        force_[0][0] += sum(fN);
        force_[1][0] += sum(fT);
        force_[2][0] += sum(fP);
        moment_[0][0] += sum(Md^fN);  // <- Line 382
        moment_[1][0] += sum(Md^fT);
        moment_[2][0] += sum(Md^fP);
    }
    else
    { … }
}

What are bins by the way? 🤔 I presume the C++ bitwise XOR operator ^ is overloaded somewhere to be used as cross product. I would be grateful for any hints where to find it in the source code so that I can have a look at how it's implemented.

I couldn't find any relevant forum posts and I built a very simple test case that is simple enough to be calculated analytically in a LibreOffice spreadsheet. If the forces are mostly integers, the forces function object outputs forces and moments that are perfectly perpendicular, with a zero scalar product and which are the same as those calculated with the spreadsheet. As soon as the forces or the face centre coordinates become floating point values, they are not perpendicular any more. Moving the CofR close to the centre of pressure, where the moments should become zero, doesn't help much.

I took one of my cases and generated the face centres in ParaView and exported the relevant data to a spreadsheet for an entire patch with about 48k faces. Apart from a general discrepancy that I found between calculating forces via function object or with ParaView, the result of the spreadsheet is in good agreement, including the total pressure force and moment not being perpendicular. Interestingly, the scalar product of any single face force and moment is zero, but once the forces and moments are summed up, the total vectors deviate from being perpendicular.

I thought it was due to a loss of accuracy of floating point numbers and so I reduced the number of faces to two. With only two faces and a 15-digit precision of LibreOffice Calc, the total forces and moment vectors are still not perpendicular. The values of forces, face centre coordinates and moments for those two faces vary by a single order of magnitude. I believe the problem lies with the summation and multiplication of a large quantity of small numbers, where the sum becoming large leads to loss of low-order significands. I wonder if there's any way to make this more accurate.

Any comments, accounts of similar experiences or suggestions are appreciated.

Thanks, Claudio

[Ship Designer]

Quick update, apparently the resultant force and moment are not perpendicular if the force system is non-planar, non-concurrent and non-parallel. It can be thus assumed that if OpenFOAM outputs such results, they are not due to numerical loss of precision but rather correctly reflect the force system as it is being created by the flow field. I have performed various calculations with LibreOffice and Octave with similar force systems which lead to the resultant force and moment not being perpendicular to one another.