[federicabi]

Hi!

I'm pretty new with openFOAM and with C++, and I'm having troubles with a customized restraint that I modified from linearSpring.
These are my .C and .H files

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "variableLinearSpring.H"
#include "addToRunTimeSelectionTable.H"
#include "sixDoFRigidBodyMotion.H"

#include "forces.H"
#include "IOobject.H"
#include "dictionary.H"
#include "Time.H"
#include "IOmanip.H"
#include "Pstream.H"

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
namespace sixDoFRigidBodyMotionRestraints
{
    defineTypeNameAndDebug(variableLinearSpring, 0);

    addToRunTimeSelectionTable
    (
        sixDoFRigidBodyMotionRestraint,
        variableLinearSpring,
        dictionary
    );
}
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::sixDoFRigidBodyMotionRestraints::variableLinearSpring::variableLinearSpring
(
    const word& name,
    const dictionary& sDoFRBMRDict   
)
:
    sixDoFRigidBodyMotionRestraint(name, sDoFRBMRDict),
    refAttachmentPt_(),
    alfaini_(),
    forzaxini_()
    
{
    read(sDoFRBMRDict);
}


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::sixDoFRigidBodyMotionRestraints::variableLinearSpring::~variableLinearSpring()
{}


// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //

void Foam::sixDoFRigidBodyMotionRestraints::variableLinearSpring::restrain
(
    const sixDoFRigidBodyMotion& motion,
    vector& restraintPosition,
    vector& restraintForce,
    vector& restraintMoment
   
) const
{
    static bool firstIt_ = true;
    scalar alfa;
    scalar forzax;
    scalar forzaxPrevIter;
    scalar alfaPrevIter;
    
  
    if (firstIt_)
    {
    alfa = (alfaini_ * 3.1415 / 180);
    forzax = forzaxini_;
    }
    else
    {
      
      FILE *forze;
      scalar timestep;
      scalar fxp;
      scalar fyp;
      scalar fzp;
      scalar fxv;
      scalar fyv;
      scalar fzv;
      scalar fxpo;
      scalar fypo;
      scalar fzpo;
      scalar mxp;
      scalar myp;
      scalar mzp;
      scalar mxv;
      scalar myv;
      scalar mzv;
      scalar mxpo;
      scalar mypo;
      scalar mzpo;
      forze = fopen("postProcessing/forces/0.3/forces.dat", "r");

      scanf("%*[^\n]");   /* Skip to the End of the Line */
    
      for (int i = 1; i <= 3; i++)
      {
      scanf("%*1[\n]");   /* Skip Three Newline */
      }

      fscanf(forze, "%le    ((%le %le %le) (%le %le %le) (%le %le %le)) ((%le %le %le) (%le %le %le) (%le %le %le))", &timestep, &fxp, &fyp, &fzp, &fxv, &fyv, &fzv, &fxpo, &fypo, &fzpo, &mxp, &myp, &mzp, &mxv, &myv, &mzv, &mxpo, &mypo, &mzpo);

      fclose(forze);
      
      forzaxPrevIter = fxp + fxv + fxpo;
      
      alfaPrevIter = (alfaini_ * 3.1415 / 180) + asin(motion.orientation()[3]);
      
      alfa = alfaPrevIter;
      forzax = forzaxPrevIter;
    }
    
    restraintPosition = motion.transform(refAttachmentPt_);

    vector v = motion.velocity(restraintPosition);

    restraintForce = vector(forzax,0,forzax*tan(alfa));

    restraintMoment = vector::zero;

    if (motion.report())
    {
        Info<< " attachmentPt" << restraintPosition
            << " force " << restraintForce
            << endl;
    }
        
    firstIt_ = false;
}

      
bool Foam::sixDoFRigidBodyMotionRestraints::variableLinearSpring::read
(
    const dictionary& sDoFRBMRDict
)
{
    sixDoFRigidBodyMotionRestraint::read(sDoFRBMRDict);

    sDoFRBMRCoeffs_.lookup("refAttachmentPt") >> refAttachmentPt_;
    sDoFRBMRCoeffs_.lookup("alfaini") >> alfaini_;
    sDoFRBMRCoeffs_.lookup("forzaxini") >> forzaxini_;
    

    return true;
}


void Foam::sixDoFRigidBodyMotionRestraints::variableLinearSpring::write
(
    Ostream& os
) const
{
    os.writeKeyword("refAttachmentPt")
        << refAttachmentPt_ << token::END_STATEMENT << nl;

    os.writeKeyword("alfaini")
        << alfaini_ << token::END_STATEMENT << nl;

    os.writeKeyword("forzaxini")
        << forzaxini_ << token::END_STATEMENT << nl;
    
    
}

// ************************************************************************* //

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

Class
    Foam::sixDoFRigidBodyMotionRestraints::variableLinearSpring

Description
    sixDoFRigidBodyMotionRestraints model. Variable Linear spring.

SourceFiles
    variableLinearSpring.C

\*---------------------------------------------------------------------------*/

#ifndef variableLinearSpring_H
#define variableLinearSpring_H

#include "sixDoFRigidBodyMotionRestraint.H"
#include "point.H"


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

namespace sixDoFRigidBodyMotionRestraints
{

/*---------------------------------------------------------------------------*\
                          Class variableLinearSpring Declaration
\*---------------------------------------------------------------------------*/

class variableLinearSpring
:
    public sixDoFRigidBodyMotionRestraint
{
    // Private data

        //- Punto di applicazione del traino (punto elica)
        point refAttachmentPt_;

        //- Angolo di inclinazione della linea d'assi
    scalar alfaini_;
    
    //- Forza di resistenza iniziale
    scalar forzaxini_;
    
    
public:

    //- Runtime type information
    TypeName("variableLinearSpring");


    // Constructors

        //- Construct from components
        variableLinearSpring
        (
            const word& name,
            const dictionary& sDoFRBMRDict
        );

        //- Construct and return a clone
        virtual autoPtr<sixDoFRigidBodyMotionRestraint> clone() const
        {
            return autoPtr<sixDoFRigidBodyMotionRestraint>
            (
                new variableLinearSpring(*this)
            );
        }


    //- Destructor
    virtual ~variableLinearSpring();


    // Member Functions

        //- Calculate the restraint position, force and moment.
        //  Global reference frame vectors.
        virtual void restrain
        (
            const sixDoFRigidBodyMotion& motion,
            vector& restraintPosition,
            vector& restraintForce,
            vector& restraintMoment
        
        ) const;
    
        
        //- Update properties from given dictionary
        virtual bool read(const dictionary& sDoFRBMRCoeff);

        //- Write
        virtual void write(Ostream&) const;
};


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace solidBodyMotionFunctions
} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //

I compiled it without errors and I applied it to a simulation with interDyMFoam solver.
I run the code and it runs for the first iteration and then it stops without any error warning, just it remains like below for ever and I don't know why.

/*---------------------------------------------------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  3.0.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
Build  : 3.0.1-d8a290b55d28
Exec   : interDyMFoam
Date   : Sep 14 2016
Time   : 14:17:37
Host   : "rd-cfd"
PID    : 45770
Case   : /home/jolly/OpenFOAM/Stage_Federica/AZ95rph/AZ95rph_snappyHexMesh/AZ95rph_snappyHexMesh_2DoFsenzaForze_26kn
nProcs : 1
sigFpe : Enabling floating point exception trapping (FOAM_SIGFPE).
fileModificationChecking : Monitoring run-time modified files using timeStampMaster
allowSystemOperations : Allowing user-supplied system call operations

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0.3

Selecting dynamicFvMesh dynamicMotionSolverFvMesh
Selecting motion solver: sixDoFRigidBodyMotion
Selecting sixDoFSolver Newmark
Translational constraint tensor (0 0 0 0 0 0 0 0 1)
Rotational constraint tensor (0 0 0 0 1 0 0 0 0)

PIMPLE: Operating solver in PISO mode

Reading field p_rgh

Reading field U

Reading/calculating face flux field phi

Reading transportProperties

Selecting incompressible transport model Newtonian
Selecting incompressible transport model Newtonian
Selecting turbulence model type RAS
Selecting RAS turbulence model kOmegaSST
Selecting patchDistMethod meshWave
kOmegaSSTCoeffs
{
    alphaK1         0.85;
    alphaK2         1;
    alphaOmega1     0.5;
    alphaOmega2     0.856;
    gamma1          0.555556;
    gamma2          0.44;
    beta1           0.075;
    beta2           0.0828;
    betaStar        0.09;
    a1              0.31;
    b1              1;
    c1              10;
    F3              false;
}


Reading g

Reading hRef
Calculating field g.h

No MRF models present

No finite volume options present

GAMG:  Solving for pcorr, Initial residual = 1, Final residual = 0.000767548, No Iterations 9
GAMG:  Solving for pcorr, Initial residual = 0.000969538, Final residual = 8.08389e-07, No Iterations 5
time step continuity errors : sum local = 4.04061e-11, global = -3.40519e-12, cumulative = -3.40519e-12
Reading/calculating face velocity Uf

Courant Number mean: 0.00117436 max: 2.90517

Starting time loop

forces forces:
    Not including porosity effects

Courant Number mean: 0.00117436 max: 2.90517
Interface Courant Number mean: 1.75985e-06 max: 1.32479
deltaT = 0.000178571
Time = 0.300179

PIMPLE: iteration 1

Restraint externalForces:  attachmentPt(1.58696 1.05 -1.73213) force (-55645 0 -2902.51)
6-DoF rigid body motion
    Centre of rotation: (9.3 0 3.18594)
    Centre of mass: (9.3 0 3.18594)
    Orientation: (1 0 -0.000808245 0 1 0 0.000808245 0 1)
    Linear velocity: (0 0 0.152943)
    Angular velocity: (0 0.0157576 0)
Execution time for mesh.update() = 9.18 s
GAMG:  Solving for pcorr, Initial residual = 1, Final residual = 0.000991738, No Iterations 11
GAMG:  Solving for pcorr, Initial residual = 0.000365091, Final residual = 7.69852e-07, No Iterations 4
time step continuity errors : sum local = 9.41351e-12, global = -6.12812e-13, cumulative = -4.01801e-12
smoothSolver:  Solving for alpha.water, Initial residual = 4.79871e-08, Final residual = 1.22406e-16, No Iterations 3
Phase-1 volume fraction = 0.545033  Min(alpha.water) = -2.91408e-11  Max(alpha.water) = 1.00003
MULES: Correcting alpha.water
MULES: Correcting alpha.water
MULES: Correcting alpha.water
Phase-1 volume fraction = 0.545033  Min(alpha.water) = -3.92276e-10  Max(alpha.water) = 1.00003
GAMG:  Solving for p_rgh, Initial residual = 0.00530249, Final residual = 4.33916e-06, No Iterations 10
GAMG:  Solving for p_rgh, Initial residual = 0.000367444, Final residual = 3.22317e-07, No Iterations 5
time step continuity errors : sum local = 4.18558e-12, global = -2.79785e-13, cumulative = -4.29779e-12
GAMG:  Solving for p_rgh, Initial residual = 9.06052e-05, Final residual = 5.46971e-08, No Iterations 6
GAMG:  Solving for p_rgh, Initial residual = 4.83781e-05, Final residual = 7.13157e-09, No Iterations 7
time step continuity errors : sum local = 9.60132e-14, global = -4.5178e-16, cumulative = -4.29824e-12
smoothSolver:  Solving for omega, Initial residual = 1.11214e-05, Final residual = 5.23818e-12, No Iterations 1
smoothSolver:  Solving for k, Initial residual = 6.08195e-08, Final residual = 5.60729e-12, No Iterations 1
ExecutionTime = 219.07 s  ClockTime = 222 s

forces forces output:
    sum of forces:
        pressure : (-1.49126e+06 -1.20981e+07 2.89098e+07)
        viscous  : (-16593.9 2035.3 -26.9414)
        porous   : (0 0 0)
    sum of moments:
        pressure : (-9.04111e+06 3.01946e+07 89945)
        viscous  : (5158.34 53168.4 42987.5)
        porous   : (0 0 0)

fieldMinMax minmaxdomain output:
    min(p) = -279303 at location (0.139124 0.399286 -1.14294)
    max(p) = 787024 at location (5.64798 0.333083 -1.02371)
    min(U) = (-41.2283 -18.1915 -41.8936) at location (18.8947 2.47892 1.3215)
    max(U) = (5.22655 7.86842 8.85684) at location (21.4319 0.100175 0.228842)

Courant Number mean: 0.000873884 max: 2.21732
Interface Courant Number mean: 1.31434e-06 max: 0.781705
deltaT = 0.000200437
Time = 0.300379

PIMPLE: iteration 1

Can anybody help me to solve my problem?
Thanks in advance

UPDATE: I tried to remove the if cycle and the firstIt dependancy and the simulation get stuck at the first iteration with no errors. I assume that the problem is about the reading of the forces.dat file but I have no idea about how to solve this problem.

Federica

[thanosniotis]

hello Federica,
Did you manage to solve this?
BW
Thanos