[aerosjc]

I use OpenFoam to simulate one particle in fluid. I used the “Spheredrag” and I would like to output the drag force on the particle. First I added some codes in the SphereDragForce.C, which is in /lagrangian/intermediate/submodels/Kinematic/ParticleForces/Drag/SphereDrag
The following is the code and the red lines are that I added by myself.

#include "SphereDragForce.H"
#include <iostream>
#include <fstream>
// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

template<class CloudType>
Foam::scalar Foam::SphereDragForce<CloudType>::CdRe(const scalar Re) const
{
    if (Re > 1000.0)
    {
        return 0.424*Re;
    }
    else
    {
        return 24.0*(1.0 + 1.0/6.0*pow(Re, 2.0/3.0));
    }
}


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

template<class CloudType>


Foam::SphereDragForce<CloudType>::SphereDragForce
(
    CloudType& owner,
    const fvMesh& mesh,
    const dictionary& dict
)
:
    ParticleForce<CloudType>(owner, mesh, dict, typeName, false)
{}


template<class CloudType>
Foam::SphereDragForce<CloudType>::SphereDragForce
(
    const SphereDragForce<CloudType>& df
)
:
    ParticleForce<CloudType>(df)
{}


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

template<class CloudType>
Foam::SphereDragForce<CloudType>::~SphereDragForce()
{}


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

ofstream outputfile1("dragforce.dat", std::ios::out );

template<class CloudType>
Foam::forceSuSp Foam::SphereDragForce<CloudType>::calcCoupled
(
    const typename CloudType::parcelType& p,
    const scalar dt,
    const scalar mass,
    const scalar Re,
    const scalar muc
) const
{
    forceSuSp value(Zero, 0.0);

    value.Sp() = mass*0.75*muc*CdRe(Re)/(p.rho()*sqr(p.d()));
	   
	outputfile1 << "drag " << value.Sp()  << endl;
	return value;
}

When I do wmake libso, the following errors appear:

Make/linux64GccDPInt32Opt/parcels/derived/basicThermoParcel/makeBasicThermoParcelSubmodels.o:(.bss+0xb80): multiple definition of `outputfile1'
Make/linux64GccDPInt32Opt/parcels/derived/basicKinematicParcel/makeBasicKinematicParcelSubmodels.o:(.bss+0xb40): first defined here
Make/linux64GccDPInt32Opt/parcels/derived/basicReactingParcel/makeBasicReactingParcelSubmodels.o:(.bss+0xd00): multiple definition of `outputfile1'
Make/linux64GccDPInt32Opt/parcels/derived/basicKinematicParcel/makeBasicKinematicParcelSubmodels.o:(.bss+0xb40): first defined here
Make/linux64GccDPInt32Opt/parcels/derived/basicReactingMultiphaseParcel/makeBasicReactingMultiphaseParcelSubmodels.o:(.bss+0xec0): multiple definition of `outputfile1'
Make/linux64GccDPInt32Opt/parcels/derived/basicKinematicParcel/makeBasicKinematicParcelSubmodels.o:(.bss+0xb40): first defined here
Make/linux64GccDPInt32Opt/parcels/derived/basicKinematicMPPICParcel/makeBasicKinematicMPPICParcelSubmodels.o:(.bss+0xdc0): multiple definition of `outputfile1'
Make/linux64GccDPInt32Opt/parcels/derived/basicKinematicParcel/makeBasicKinematicParcelSubmodels.o:(.bss+0xb40): first defined here
collect2: error: ld returned 1 exit status
make: *** […/platforms/linux64GccDPInt32Opt/lib/libZHOnelagrangianIntermediate.so] Error 1

So I think maybe we cannot add the codes in SphereDragForce.C to output the drag force on particle. Then I added some codes (red color) in the creatFields.H and here are the codes:

#include "readGravitationalAcceleration.H"

word continuousPhaseName
(
    IOdictionary
    (
        IOobject
        (
            "transportProperties",
            runTime.constant(),
            mesh,
            IOobject::MUST_READ
        )
    ).lookup("continuousPhaseName")
);

Info<< "Reading field U\n" << endl;
volVectorField Uc
(
    IOobject
    (
        IOobject::groupName("U", continuousPhaseName),
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

Info<< "Reading field p\n" << endl;
volScalarField p
(
    IOobject
    (
        "p",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);


Info<< "Reading/calculating continuous-phase face flux field phic\n"
    << endl;

surfaceScalarField phic
(
    IOobject
    (
        IOobject::groupName("phi", continuousPhaseName),
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        IOobject::AUTO_WRITE
    ),
    linearInterpolate(Uc) & mesh.Sf()
);

label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, pimple.dict(), pRefCell, pRefValue);
mesh.setFluxRequired(p.name());

Info<< "Creating turbulence model\n" << endl;

singlePhaseTransportModel continuousPhaseTransport(Uc, phic);

dimensionedScalar rhocValue
(
    IOobject::groupName("rho", continuousPhaseName),
    dimDensity,
    continuousPhaseTransport.lookup
    (
        IOobject::groupName("rho", continuousPhaseName)
    )
);

volScalarField rhoc
(
    IOobject
    (
        rhocValue.name(),
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    rhocValue
);

volScalarField muc
(
    IOobject
    (
        IOobject::groupName("mu", continuousPhaseName),
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    rhoc*continuousPhaseTransport.nu()
);

Info << "Creating field alphac\n" << endl;
// alphac must be constructed before the cloud
// so that the drag-models can find it
volScalarField alphac
(
    IOobject
    (
        IOobject::groupName("alpha", continuousPhaseName),
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedScalar("0", dimless, 0)
);

word kinematicCloudName("kinematicCloud");
args.optionReadIfPresent("cloudName", kinematicCloudName);


volScalarField spheredrag
(
IOobject
(
"SphereDragForce",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
);



Info<< "Constructing kinematicCloud " << kinematicCloudName << endl;
basicKinematicTypeCloud kinematicCloud
(
    kinematicCloudName,
    rhoc,
    Uc,
    muc,
    g
);

// Particle fraction upper limit
scalar alphacMin
(
    1.0
  - readScalar
    (
        kinematicCloud.particleProperties().subDict("constantProperties")
       .lookup("alphaMax")
    )
);

// Update alphac from the particle locations
alphac = max(1.0 - kinematicCloud.theta(), alphacMin);
alphac.correctBoundaryConditions();

surfaceScalarField alphacf("alphacf", fvc::interpolate(alphac));
surfaceScalarField alphaPhic("alphaPhic", alphacf*phic);

autoPtr<PhaseIncompressibleTurbulenceModel<singlePhaseTransportModel>>
continuousPhaseTurbulence
(
    PhaseIncompressibleTurbulenceModel<singlePhaseTransportModel>::New
    (
        alphac,
        Uc,
        alphaPhic,
        phic,
        continuousPhaseTransport
    )
);

When I run wmake, it compile successfully.
But when simulate a case, it shows the error
cannot find file "…/0/SphereDragForce"

From function virtual Foam::autoPtr<Foam::ISstream> Foam::fileOperations::uncollatedFileOperation::rea dStream(Foam::regIOobject&, const Foam::fileName&, const Foam::word&, bool) const in file
global/fileOperations/uncollatedFileOperation/uncollatedFileOperation.C at line 505.

FOAM exiting


Thanks

[saiguruprasad]

Did you find out how to print the drag force on the particle?

[Shah Akib Sarwar]

Did any of you find a solution to this? I am trying to do the same thing, but cannot figure out how.