[upuli]

Hi



I want to move a scalar which has value 1 to be moved along x axis at a constant speed. It has no diffusion or source terms. When I use +fvm::div term to calculate the movement due to velocity it gives fractions to the results. Can some one help me on how to do it?


regards


Upuli

[upuli]

Hi


I write the conservation equation for scalar solid as below.

fvScalarMatrix solidEqn
        (
            fvm::ddt(solid)
       fvc::div(phigrate,solid)
        );

    solve(solidEqn==Sc);

When I execute the solver the value of the solid is shown with the decimals (something like 1.24 etc) due to the divergence term. But when the scalar solid is moved at a constant velocity I want to have 1 or 0 as the value of the solid . Can someone please help me to fix the issue.


Thanks

[msaravia]

I think that probably you have numerical diffusion, check the discretization scheme for the convection term. Also, the convection term should be evaluated implicitly. See scalarTransporFoam solver.

[Tobi]

Can you explain us, why you do not solve the scalar as:

solidEqn.solve();

What is 'Sc' ? In addition to the previous reply. What are your schemes and your solution tactic? Is everything converged? If you use 'Upwind' you should be between 0-1. All other schemes can lead to non-physical results. See my website and the numerical schemes topic.

Good luck.

[upuli]

Thank you very much for the reply

Sc is a zero valued scalar .Actually I forgot that it can be written as

solidEqn.solve();

When the Gauss upwind scheme is used the the solid does not seem to move. This transport equation was added to a my own solver. It gives printstack error after about 30 seconds and stops.


Thanks

[Tobi]

Without any information, you will not get any support. It is like: My cooked meal is not good, so tell my why?

If the scalar is not moving with Upwind it is an indication that something went wrong in your development. In addition, relaxing the scalar or matrix can help.

Good luck.

[upuli]

Hi


sorry for inadequate data. Now I removed all the other conservation equations from the solver and only the solid conservation equation is present in the solver. It is the icoFoam solver. The solid movement appears only when

+fvc::div(phigrate,solid)

term is used. Fv scheme used is Upwind. The solver specifications are as below

    solid
    {
        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-05;
        relTol          0;
    }

the boundaries for velocity ugrate is 0.001 m/s in the positive x direction for all the patches.
Now the solver is running and it is able to produce 1 and 0 only.

When

+fvm::div(phigrate,solid)

term is used the solids does not move(for any of the fvschemes).

[Ben UWIHANGANYE]

Hello Foamers,

I would like to ask how the modified setFieldsDict is used in FSI simulation.

Indeed, I want to simulate flow over a cylinder(stationary and oscillating). I have used toposetDict to map that cylinder in fluid Cartesian domain.

Let: eta=1 solid region, eta=0 fluid region.
Then i have a c++ code that gives me a value of eta each time step. The eta's value alternate from 0 to 1 and from 1 to 0. C++ code prints eta values in .dat format.

I am asking if it is possible to use setFieldsDict to move (oscillating motion) that cylinder as eta changes.

What should be my code looking like? can advise how i can rotate that solid region with that eta alternation?

I will be happy to hear from you!

Regard!

[upuli]

Hi

What I used is setFields not the toposet Fields. I think both of them are similar. I modeled the translation movement. So the velocity component in x direction only. With the conditions given in previous post , my solver works correctly and it only gives 0 and 1. I think with correct representation of oscillating velocity (momentum) the movement of cylinder can be done.

I think rotation movement can be done if you know the angular velocity and taking the tangential velocity components in x and y directions.

regards

Upuli

[Ben UWIHANGANYE]

Originally Posted by upuli

Hi

What I used is setFields not the toposet Fields. I think both of them are similar. I modeled the translation movement. So the velocity component in x direction only. With the conditions given in previous post , my solver works correctly and it only gives 0 and 1. I think with correct representation of oscillating velocity (momentum) the movement of cylinder can be done.

I think rotation movement can be done if you know the angular velocity and taking the tangential velocity components in x and y directions.

regards

Upuli

Thank you for your reply,

I have attached a c++ program on this post. It is for undermped free vibration. the solid is modeled to vibrate freely. Can you help to know it can work in openFoam to rotate that cylinder.

May you share the code you used(benclaude2011@gmail.com) in your case.

Thank you!

[upuli]

HI


I modified the icoFoam solver.below is what I have used

#include "fvCFD.H"
#include "pisoControl.H"

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

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"

    pisoControl piso(mesh);

    #include "createFields.H"
    #include "initContinuityErrs.H"

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

    Info<< "\nStarting time loop\n" << endl;

    while (runTime.loop())
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        #include "CourantNo.H"

        // Momentum predictor
fvScalarMatrix solidEqn
        (
            fvm::ddt(solid)
       +fvc::div(phigrate,solid)
        );
solidEqn.solve();
        fvVectorMatrix UEqn
        (
            fvm::ddt(U)
          + fvm::div(phi, U)
          - fvm::laplacian(nu, U)
        );

        if (piso.momentumPredictor())
        {
            solve(UEqn == -fvc::grad(p));
        }

        // --- PISO loop
        while (piso.correct())
        {
            volScalarField rAU(1.0/UEqn.A());
            volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
            surfaceScalarField phiHbyA
            (
                "phiHbyA",
                fvc::flux(HbyA)
              + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
            );

            adjustPhi(phiHbyA, U, p);

            // Update the pressure BCs to ensure flux consistency
            constrainPressure(p, U, phiHbyA, rAU);

            // Non-orthogonal pressure corrector loop
            while (piso.correctNonOrthogonal())
            {
                // Pressure corrector

                fvScalarMatrix pEqn
                (
                    fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
                );

                pEqn.setReference(pRefCell, pRefValue);

                pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));

                if (piso.finalNonOrthogonalIter())
                {
                    phi = phiHbyA - pEqn.flux();
                }
            }

            #include "continuityErrs.H"

            U = HbyA - rAU*fvc::grad(p);
            U.correctBoundaryConditions();
        }

        runTime.write();

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }

    Info<< "End\n" << endl;

    return 0;
}


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

The below code may help you. But it may not be a perfect way to do it.

fvScalarMatrix etaEqn
        (
            fvm::ddt(eta)
       +fvc::div(phieta,eta)
        );
etaEqn.solve();

I think you have to write conservation for the momentum change of eta

fvVectorMatrix UetaEqn
        (
            fvm::ddt(Ueta)
          + fvm::div(pheta, Ueta)
         
        );
 solve(UetaEqn ==−ω^2*x) ;

The above is only a suggestion. May not be a perfect method to solve your problem.



p { margin-bottom: 0.1in; direction: ltr; line-height: 120%; text-align: left; }a:link { }

[Ben UWIHANGANYE]

Thank you,

I will try it and will let you know the progress.
I am asking:
no changes needed in creatField.H file? after compilation, what about in O folder fvscheme and fvsolution ?

what will be the impact of solve(UetaEqn ==−ω^2*x) in my case, will i be still in need of c++ program.


Regard!

[upuli]

Hi


You have to declare the scalar field eta and vector field Ueta in creatField.H file.


regards

[Ben UWIHANGANYE]

Dear Upuli,

I have tried to go as you advised but after compiling i have the following errors:

Making dependency list for source file myicoFoam.C
g++ -std=c++0x -m64 -Dlinux64 -DWM_ARCH_OPTION=64 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -O3 -DNoRepository -ftemplate-depth-100 -I/opt/openfoam4/src/finiteVolume/lnInclude -I/opt/openfoam4/src/meshTools/lnInclude -IlnInclude -I. -I/opt/openfoam4/src/OpenFOAM/lnInclude -I/opt/openfoam4/src/OSspecific/POSIX/lnInclude -fPIC -c myicoFoam.C -o Make/linux64GccDPInt32Opt/myicoFoam.o
myicoFoam.C:69:2: error: stray ‘\342’ in program
solve(UetaEqn ==−ω^2*x) ;
^
myicoFoam.C:69:2: error: stray ‘\210’ in program
myicoFoam.C:69:2: error: stray ‘\222’ in program
myicoFoam.C:69:2: error: stray ‘\317’ in program
myicoFoam.C:69:2: error: stray ‘\211’ in program
myicoFoam.C: In function ‘int main(int, char**)’:
myicoFoam.C:66:22: error: ‘phieta’ was not declared in this scope
+ fvm::div(phieta, Ueta)
^
myicoFoam.C:69:23: error: expected primary-expression before ‘^’ token
solve(UetaEqn ==−ω^2*x) ;
^
myicoFoam.C:69:26: error: ‘x’ was not declared in this scope
solve(UetaEqn ==−ω^2*x) ;
^
/opt/openfoam4/wmake/rules/General/transform:8: recipe for target 'Make/linux64GccDPInt32Opt/myicoFoam.o' failed
make: *** [Make/linux64GccDPInt32Opt/myicoFoam.o] Error 1



The following are the change made in myicoFoam.C:

#include "fvCFD.H"
#include "pisoControl.H"

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

int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"

pisoControl piso(mesh);

#include "createFields.H"
#include "initContinuityErrs.H"

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

Info<< "\nStarting time loop\n" << endl;

while (runTime.loop())
{
Info<< "Time = " << runTime.timeName() << nl << endl;

#include "CourantNo.H"

// Momentum predictor


fvVectorMatrix UetaEqn
(
fvm::ddt(Ueta)
+ fvm::div(phieta, Ueta)

);
solve(UetaEqn ==−ω^2*x) ;


fvScalarMatrix etaEqn
(
fvm::ddt(eta)
+fvc::div(phieta,eta)
);
etaEqn.solve();


fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
- fvm::laplacian(nu, U)
);

if (piso.momentumPredictor())
{
solve(UEqn == -fvc::grad(p));
}

// --- PISO loop
while (piso.correct())
{
volScalarField rAU(1.0/UEqn.A());
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
surfaceScalarField phiHbyA
(
"phiHbyA",
fvc::flux(HbyA)
+ fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
);

adjustPhi(phiHbyA, U, p);

// Update the pressure BCs to ensure flux consistency
constrainPressure(p, U, phiHbyA, rAU);


// Non-orthogonal pressure corrector loop
while (piso.correctNonOrthogonal())
{
// Pressure corrector

fvScalarMatrix pEqn
(
fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
);

pEqn.setReference(pRefCell, pRefValue);

pEqn.solve(mesh.solver(p.select(piso.finalInnerIte r())));

if (piso.finalNonOrthogonalIter())
{
phi = phiHbyA - pEqn.flux();
}
}

#include "continuityErrs.H"

U = HbyA - rAU*fvc::grad(p);
U.correctBoundaryConditions();

}

runTime.write();

Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
<< " ClockTime = " << runTime.elapsedClockTime() << " s"
<< nl << endl;
}

Info<< "End\n" << endl;

return 0;
}


my creatField.H:

Info<< "Reading transportProperties\n" << endl;

IOdictionary transportProperties
(
IOobject
(
"transportProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);

dimensionedScalar nu
(
"nu",
dimViscosity,
transportProperties.lookup("nu")
);

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


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

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

#include "createPhi.H"


label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
mesh.setFluxRequired(p.name());

May you help to know what is wrong, and if you can, you may let me know how your 0 folder was looking like to move that solid scalarField.


Thank you!

[upuli]

Hi


You have to declare phieta and x , since these fields are created by you ,you have to declare them. They are not originally included in OpenFOAM. When you declare x ,also you have to take the displacement vector. And sorry for the wrong representation of mathematical operators. In OpenFOAM (^operator doesnot work) Please refer the programmers guide.


It should be

fvVectorMatrix UetaEqn
        (
            fvm::ddt(Ueta)
          + fvm::div(phieta, Ueta)
         
        );
 solve(UetaEqn ==−pow(ω,2)*x) ;

[Ben UWIHANGANYE]

The remaining errors are:

g++ -std=c++0x -m64 -Dlinux64 -DWM_ARCH_OPTION=64 -DWM_DP -DWM_LABEL_SIZE=32 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -O3 -DNoRepository -ftemplate-depth-100 -I/opt/openfoam4/src/finiteVolume/lnInclude -I/opt/openfoam4/src/meshTools/lnInclude -IlnInclude -I. -I/opt/openfoam4/src/OpenFOAM/lnInclude -I/opt/openfoam4/src/OSspecific/POSIX/lnInclude -fPIC -c myicoFoam.C -o Make/linux64GccDPInt32Opt/myicoFoam.o
myicoFoam.C:69:2: error: stray ‘\342’ in program
solve(UetaEqn ==−ω*ω*x);
^
myicoFoam.C:69:2: error: stray ‘\210’ in program
myicoFoam.C:69:2: error: stray ‘\222’ in program
myicoFoam.C:69:2: error: stray ‘\317’ in program
myicoFoam.C:69:2: error: stray ‘\211’ in program
myicoFoam.C:69:2: error: stray ‘\317’ in program
myicoFoam.C:69:2: error: stray ‘\211’ in program
myicoFoam.C: In function ‘int main(int, char**)’:
myicoFoam.C:69:27: error: ‘x’ was not declared in this scope
solve(UetaEqn ==−ω*ω*x);
^
/opt/openfoam4/wmake/rules/General/transform:8: recipe for target 'Make/linux64GccDPInt32Opt/myicoFoam.o' failed
make: *** [Make/linux64GccDPInt32Opt/myicoFoam.o] Error 1


I am asking if all errors related to the displacement vector x which is not yet declared.

Can you help to know how I can declare it.


Regard!

[upuli]

Hi

You have to declare it in create fields file. Since x is the displacement in either x,y or z directions you have to refer on how to write displacement in OpenFoam with reference to a certain direction.

Good luck