Maybe minus should be a plus in FSI/Force.C

[moaxm]

https://github.com/precice/openfoam-adapter/blob/8f0ff7d27cd3a1de6a393465ef7516eb8b897794/FSI/Force.C#L59

Maybe this gives you strange tangential forces at the moment.

[derekrisseeuw]

Hi @moaxm

Thank you for the feedback. I am curious how you found this error, could you please explain your simulation setup?

[moaxm]

Hi @derekrisseeuw We have created a similar socket-based communication interface between OpenFOAM and Abaqus which we use within our coupling environment ifls (https://www.tu-braunschweig.de/ifl/simulationswerkzeuge/ifls) and in our interface we do not have a minus there.

[derekrisseeuw]

Interesting remark. As far as I can see there is no problem with the implementation in the openfoam adapter. We used that same convention as is used in the forces class in openfoam, see for reference:

force.C

Also, I checked quickly to see whether the velocity and force vector on the boundary were headed in the same direction, and they are.

Could it be that the minus sign is located somewhere else in your code? Maybe it is migrated to the part where the viscous and pressure forces are added. (or subtracted)

[moaxm]

Hi, we calculate the total force as follows:

const surfaceVectorField::Boundary& Sfb =
    mesh_.Sf().boundaryField();

tmp<volSymmTensorField> tdevRhoReff = -rho()*nu*dev(twoSymm(fvc::grad(U)));

const volSymmTensorField::Boundary& devRhoReffb
    = tdevRhoReff().boundaryField();

vectorField fN
(
    rho(p)*Sfb[patchID_]*p.boundaryField()[patchID_]
);

vectorField fT(Sfb[patchID_] & devRhoReffb[patchID_]);

The total force is simply fN+fT. The total force is then directly mapped onto the structural mesh. The relative deviation of the results with 203k fluid cells for the deflection of the flag are below 5% comparing to the reference value given in the TH paper.

With the minus we got oscillations of the flag in the FSI3-Benchmark in x-direction.

I verified the viscous force by taking a look at the fluid velocity gradient at the surface of the flag and then compared the expected direction of the force vector with the Abaqus results.

[davidscn]

The direction of the force vector in the fluid solver depends also on the implementation in the structure solver. This might be different in this case. I simulated a coupled FSI benchmark using deal.II and the current calculation for the forces. However, the direction mapping for the x-direction was consistent in this setup

[derekrisseeuw]

You've got me a bit puzzled here @moaxm. As far as I can see, our expressions are equivalent. If I am not mistaking, you are talking about the minus sign for the viscous forces, which you include yourself too, before the rho():

tmp<volSymmTensorField> tdevRhoReff = -rho()*nu*dev(twoSymm(fvc::grad(U)));

Looking at our implementation, the following lines should perform the same:

Foam::tmp<Foam::volSymmTensorField

 preciceAdapter::FSI::Force::devRhoReff(dimensionedScalar rho) const
{ ...
    return -rho * nu * dev(twoSymm(fvc::grad(U)));
}

tmp<volSymmTensorField> tdevRhoReff = devRhoReff(rho);
const volSymmTensorField::Boundary& devRhoReffb = tdevRhoReff().boundaryField();

const volScalarField& p = mesh_.lookupObject<volScalarField>("p");

Then, in the loop:

        // TODO: Extend to cover also compressible solvers
        Force_->boundaryFieldRef()[patchID] =
            Sfb[patchID] * p.boundaryField()[patchID]* rho.value();

        // Viscous forces
        Force_->boundaryFieldRef()[patchID] +=
            Sfb[patchID] & devRhoReffb[patchID];

[moaxm]

Nevermind, I initially compared with the formulation for turbulent flows where the minus is already included in a separate devReff() function, thus the misconception that there is not a minus in our case. But you are perfectly right - the laminar formulations match in both cases. Sorry for bringing this up!

[derekrisseeuw]

Thanks for your honest feedback! Good to see that some experienced FSI modellers are taking a look at the code

[MakisH]

Thank you very much, @moaxm for looking into this and @derekrisseeuw and @DavidSCN for answering. We looked into this very quickly with @moaxm and I asked him to open an issue, but we did not have the time to properly test it. Eitherways, it's really great to know which parts of the code are correct!

Please, everybody, if you see anything else that looks out-of-place, don't hesitate and let us know!