CFD Online Discussion Forums - OpenFOAM Programming and Development

MULES::limiter

santiagomarquezd — Wed, 08 Sep 2010 01:41:46 GMT

Hi all, I've been studying MULES for a while and can understand the main ideas of the integrator, but I have some problems interpreting the idea behind the limiter. Usually one uses limiters to limit the flux before integrating temporally. In MULES a flux is given to the integrator, for example, MULES::explicitSolve, this flux is already limited by the advection schemes used to assemble the flux.
Why another limiter function is applied to the flux before integrating temporally?
Is it to complete a work partially done by interfaceCompression scheme?
Is there some reference, at least similar, to limiter used in MULES? Is it loosely based on what???

All ideas are useful. Regards.

For loops in OpenFoam

balkrishna — Mon, 06 Sep 2010 11:44:14 GMT

I am generalising the twoPhase Euler Foam to handle n Phases . In the file pEqn.H I want to sum up the phase fraction of all the dispersed phases to get the dispersed phase . To do so I programmed the following :

Code:

 PtrListalphaf(nPhases);                                                                                               

                                                                                                                                              

    //just an initialization                                                                                                                  

    surfaceScalarField af ;                                                                                                                   

    surfaceScalarField betaf = scalar(1) - af;                                                                                                

                                                                                                                                              

    forAll(alphaf,aPh)                                                                                                                        

    {                                                                                                                                         

        alphaf.set(aPh,fvc::interpolate(phases[aPh].alpha()));                                                                                

        af= af + alphaf[aPh];                                                                                                                 

    }                                                                                                                                         

    betaf = scalar(1) - af ;

On compilation I got the following error :

Code:

error: no matching function for call to �Foam::GeometricField::GeometricField()�

/home/ifmg/OpenFOAM/OpenFOAM-1.7.0/src/OpenFOAM/lnInclude/GeometricField.C:652: note: candidates are: Foam::GeometricField::GeometricField(const Foam::IOobject&, const Foam::GeometricField&, const Foam::wordList&) [with Type = double, PatchField = Foam::fvsPatchField, GeoMesh = Foam::surfaceMesh]

scalar + volScalarField

akidess — Fri, 03 Sep 2010 15:11:39 GMT

Hello everyone,

I'm getting a Floating Point Exception on the following code fragment:
Foam::log(scalar(1) + (a_O*Kseg_O ))

a_O is a dimensionedScalar with value 0, Kseg_O is a volScalarField. The value in the brackets should never be zero thanks to scalar(1), yet still I get a floating point exception. Does the code not add scalar(1) to each element of the volScalarField?

Thanks,

Anton

Interpolation from volume field onto a surface

dasch — Thu, 02 Sep 2010 16:44:02 GMT

Hi all,

I'm trying to interpolate the values from a volumeScalarField onto a certain surface patch (this is going to be a coupling between OF and another solver - at some point). Is there a boundary condition for this kind of problem? If not, what do I have to do? I'm really stuck on this.

Thanks a lot for your help
Daniela

How reactingParcel updates it date after reaction?

Edison_Ge — Thu, 02 Sep 2010 12:50:03 GMT

Hi, guys, I'm crying out desperately for help. I'm developing a new reacting particle class. But how can I update its particle compositions by chemical reaction?
I think in the existing reactingParcel class function calc() is used to do it, but I didn't see the correspondent reactingCloud call this calc() function. Can someone explain a little about how reactingParcel reactes?
Edison

Adding acceleration term in icoDyMFoam

rieuk — Tue, 31 Aug 2010 14:51:01 GMT

Hey guys, real simple question. My moving reference frame is accelerating -0.5*pi*A*T*sin(2*pi*t/T) in x-dirn and -sqrt(3)/2*pi*A*T*sin(2*pi*t/T) in y-dirn. A, T, pi are constants, t is time. I'm trying to add these to the momentum equation but it doesn't compile. Here's what I've got:

Ueqn.H

Code:

volScalarField pi = 3.1415926535897932384626433832795028841971693993751;

volVectorField acc = -pi*0.1091911*sin(2*pi*runtime.time().value()/1.091911)*(0.5*vector(1,0,0)+pow(3,0.5)*vector(0,1,0));

    

fvVectorMatrix UEqn

    (

        fvm::ddt(U)

      + fvm::div(phi, U)

      - fvm::laplacian(nu, U)

      - acc

    );



    if (momentumPredictor)

    {

        solve(UEqn == -fvc::grad(p));

    }

What's wrong with it??

scale(a, b) function

ozgur — Tue, 31 Aug 2010 13:05:27 GMT

Hi,

In the Programmer�s Guide there is a description of a function named as scale:

Quote:

As the name suggests, the scale function is a tool for scaling the components of one tensor by the components of another tensor of the same rank. It is evaluated as the product of corresponding components of 2 tensors, e.g., scaling vector a by vector b would produce vector c whose components are ci = scale(a, b) = (a1b1, a2b2, a3b3)

I have tried to use it however, I get the following error during compilation:

Code:

error: �scale� was not declared in this scope

My code is simple:

Code:

#include "fvCFD.H"

...

vector a(1,2,3), b(4,5,6), c;

c=scale(a,b)

What can I do about this?

Thanks,

Ozgur

If-else condition over the whole mesh

vkrastev — Tue, 31 Aug 2010 10:33:45 GMT

Hi everybody, I have some troubles with finding the correct way to modify a source .C file in order to seek for positive/negative scalar values over the whole range of mesh cells. Namely, I'm trying to implement an improved version of the k-omega turbulence model, in which there is a correction of some of the closure coefficients that implies an if-else "check" over a field of scalar values. The following is the part of the source code containing the if-else condition (in the form I trivially suppose it could be correct):

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

tmp kOmegaSI::kiKappa() const
{
return (1/(pow(omega_,3.0)))*((fvc::grad(k_))&(fvc::grad(ome ga_)));
}

tmp kOmegaSI::fBetaStar() const
{
{
forAll(mesh.cells(),celli)
if (kiKappa()[celli] > scalar(0.0)) {
return ((scalar(680.0)+sqr(kiKappa()[celli]))/(scalar(400.0)+sqr(kiKappa()[celli])));
}
else {
return scalar(1.0);
}
}
}

And here it is the error message which appears after running the wmake utility:

kOmegaSI_1/kOmegaSI.C:56: error: �mesh� was not declared in this scope
kOmegaSI_1/kOmegaSI.C:57: error: no match for �operator[]� in �Foam::incompressible::RASModels::kOmegaSI::kiKapp a() const()[celli]�
kOmegaSI_1/kOmegaSI.C:58: error: no match for �operator[]� in �Foam::incompressible::RASModels::kOmegaSI::kiKapp a() const()[celli]�
kOmegaSI_1/kOmegaSI.C:58: error: no match for �operator[]� in �Foam::incompressible::RASModels::kOmegaSI::kiKapp a() const()[celli]�
kOmegaSI_1/kOmegaSI.C:61: error: conversion from �double� to non-scalar type �Foam::tmp >� requested

So, can anybody give me some hints about my problem? Of course, if necessary, I could add additional information about the .C (and the corresponding .H header) file.

Thank you in advance

Adding U*U as a source term

ozgur — Mon, 30 Aug 2010 20:45:30 GMT

Hi Foamers,

I am trying to solve the following equation:

dU/dt = grad(p) - 0.5*(f/L)*(U*U)

here, U(velocity vector), p (pressure), and L is length, (f is dimensionless scalar).

In fact, I made an attempt and it seems working for 1D .. However, I can not be very sure about the approach that I use (especially for 2D and 3D problems), and I need some comments from more experienced people.

I constructed my equations as follows:

Code:

        ...



        vector identity(1,1,1);



        ...

    

        volVectorField gradp = fvc::grad(p);



        ...



        fvVectorMatrix UEqn

        (

                 fvm::ddt(U) + fvm::Sp((U & identity)/2*(1/L+f/D), U)

        );



        ...



        solve(UEqn == -fvc::grad(p));



        ...

The reason that I am not sure is about the term "U*U" in the original equation. Since U is a vector, U*U will produce a matrix, and therefore it was not possible to write "U*U" as it is to the code as a source term.

Instead, assuming that only the i component of the vectors will be used while solving Ux, I wrote "U & identity" to obtain the value of the x component of U and put it in to the first parameter place of Sp function.

This works for 1D but I cannot be sure about my assumption about the behavior of U & identity expression in 2D and 3D problems..

In short, say that U(1, 2, 3) then will "U & identity" produce 1+2+3=6 in ALL Ux, Uy and Uz solutions OR will it produce 1 in Ux solution, 2 in Uy solution and 3 in Uz solution??

Kind regards,

Ozgur

Membrane and Darcy's law

NewFoamer — Sun, 29 Aug 2010 11:41:40 GMT

Hello,

I'm new to OpenFOAM, and as a starting project, I'm trying to simulate the crossflow over a permeable membrane. The setup is roughly sketched in this figure:

For the flow through the membrane, I simply wish to put an "internal" boundary condition on the membrane, setting:

velocity_through_membrane=k*Delta_P_over_membrane

So, my question is; how do I implement such an internal boundary condition?

Also, what would be a proper boundary condition for the pressure at the membrane?

Adding a term that includes time

ozgur — Fri, 27 Aug 2010 14:05:31 GMT

Hi,

A very simple question: In OpenFOAM, how can I add a term that includes time into my equations?

My first attempt which doesn't work is as follows:

solve
(
fvm::ddt(T)==-exp(-time())
);

Thanks for any comments,

Ozgur

Particle position each timestep (during simulation)

lindner — Fri, 27 Aug 2010 10:32:20 GMT

Hi, i have the following problem:
I need to know how many particles are there in the zone 1 and zone 2 of a room.

An example:
10 meter long room (x direction), while the first 5 meters are the zone 1 and the other 5 meters are zone 2.
The particle position file outputs the coordinate (x y z) of each particle, so it will be like:
if x < 5, sum to zone 1,
else, sum to zone 2

The actual output while solving the uncoupledKinematicParcelFoam is

Cloud: kinematicCloud
Total number of parcels added = 1210
Total mass introduced = 1.01368723e-17
Current number of parcels = 1210
Current mass in system = 1.01368723e-17
ExecutionTime = 0.61 s ClockTime = 1 s

and my objective is to add the following
Current number of parcels in Zone 1 = xxx
Current number of parcels in Zone 2 = yyy

Is it possible to make it in the solver?
Or it's better to make a separate code to extract this information?
By the way, i don't know c++ programming =(

Thank you, and if anyone needs I can send the room from Annex 20 with the uncoupledKinematicParcelFoam working.

PS: Using OpenFOAM 1.7.x

OpenFOAM MPI coupling

Ford Prefect — Mon, 23 Aug 2010 12:00:37 GMT

Hello,

I am not a current user of OpenFOAM but I might be soon ;)

I have an external software that I wish to call at certain times for some of the computational cells. I wish to distribute the calls to this external software evenly in an array using MPI. I wonder if it is possible to invoke MPI commands in user defined functions (or whatever it is called in OF)?

Best Regards

Reference pressure

OFU — Thu, 19 Aug 2010 08:24:07 GMT

Hello,

for my computations I need to set a reference pressure in one cell of the outer wall (named oustide) of my geometry. But how can I find out which cells are part of the outer wall und how they are numbered, so that I can choose one from these cells? I had a look in the five files boundary, faces, neighbour, owner, points, but that didn't help me.
And how can I set the reference value in one cell? By using pRefCell and pRefValue?

Would be nice, if anyone could help me.

Thanks,
Jennifer

Uncoupled Eulerian particle transport

peterdhicks — Tue, 17 Aug 2010 15:30:14 GMT

Hello

I'm trying to calculate the evolution of particles within a steady pre-calculated air flow using an Eulerian model. For the problem I'm trying to solve (droplet transport through air towards impact with a solid body in aircraft icing), the volume of droplets is very low, so there's no coupling to the air, and no coupling between droplets (so droplet momentum conservation involves a drag term, but no pressure).

Firstly, can this case be done with twoPhaseEulerFoam? From my limited understanding of the source code this seems to require both phases to be unsteady.

I've been trying to extend scalarTransportFoam to add a droplet momentum conservation equation and I haven't yet got satisfactory results. The guts of my solver are

Code:

            solve

            (

                fvm::ddt(alpha)

              + fvm::div(phiDrop, alpha)

            );



            volVectorField Ur = U - UDrop;

            volScalarField magUr = mag(Ur);



            fvVectorMatrix UDropEqn

            (

                fvm::ddt(UDrop)

              + fvm::div(phiDrop, UDrop)

              - 0.03*magUr*Ur/diamDrop

            );



            solve (UDropEqn);

and I was wondering if I was going about it the right way. The drag model is initial very basic for simplicity.

One of the problems I don't seem to be able to get right are the boundary conditions for the droplet velocity (UDrop) and volume fraction (alpha) on the solid body. [1] Gives a boundary condition -alpha*(UDrop dot n) > 0 on the solid surface, but I'm not sure how to replicate this in OpenFOAM. I've been trying inletOutlet and outletInlet conditions, so far without success.

Any thoughts/comments/assistance would be greatly appreciated.

Peter

[1] http://dspace.uta.edu/bitstream/hand...pdf?sequence=1