[mdjames]

Hoping someone has been down this dark and winding trail before. I've seen some modifications to other codes, for instance icoFoam and channelFoam but not this guy...

I'm a Fortran-convert new FOAM user and am having a good amount of trouble, trying to add a source term to the pimpleFoam solver.

I take it it's not as easy as adding a constant to:

tmp<fvVectorMatrix> UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
+ turbulence->divDevReff(U)
);

in the UEqn.H file at $FOAM_SOLVERS/incompressible/pimpleFoam

I'd like to have a driven half-channel flow with cyclic b.c's all around. I'm trying to get away from using inlet/outlet patches. I was using fixedValue for pressure and pressureInletVelocity for U but wasn't getting the behavior I need. I'm looking for a little less assumption on the code's part (hence the desire for the body force)

Thanks a TON for any help!

[alberto]

In principle, you can add the forcing term (it is a constant vector, not a scalar constant, since you are formally imposing a constant mean pressure gradient) to the UEqn. There is nothing that prevents you from doing that.

If you want to do it in a more careful way, you might want to include it in the flux instead than in the UEqn, before solving for the pressure equation. The treatment is identical to the one used for the gravity, for example, in bubbleFoam, and it requires some more modifications to the code.

Best,

[mdjames]

Ah, great; she's off and running. Thanks for the quick reply!

I created an H file (class?) to read a dict file which specifies the body force as a vector and its dimensions. After that, I stuck the force into the UEqn.

If you'll pardon my ignorance, what do you mean by including it in the flux?

[mdjames]

Oh, and why is that preferable to including the term in the UEqn?

[alberto]

Including force terms in the flux is a way to improve the robustness of the algorithm in some cases (for example when the force is related to the density, like for the gravity).

If you write your momentum equation as

ddt(U) + div(UU) = div(tau) - grad(p) + F

where F is your forcing term, and then you write in semi-discrete form you get

A*U = H + F

where I put everything I do not want to carry around in H. Now, the predicted velocity is

U = H/A + F/A

and the flux, obtained interpolating this on the cell faces (indicated by _f) is

phi = (H/A)_f . S + (F/A)_f . S

In your case F is constant, so

phi = (H/A)_f . S + F(1/A)_f . S

In OpenFOAM notation, this means that you should add

surfaceScalarField rUAf = fvc::interpolate(rUAf)

surfaceScalarField phiForcing = rUAf*(F & mesh.Sf())

to your flux, after adjustPhi. At this point however you have to modify the correction step for the velocity too (the flux is fine), consistently with the new velocity predictor:

U += fvc::reconstruct(phiForcing) - rUA*fvc::grad(p);

On a side note, you probably won't need this in your case, but it is better to know

Hopefully I did not introduce errors while typing

Best,

[alberto]

Originally Posted by mdjames

Oh, and why is that preferable to including the term in the UEqn?

Numerically this is done to improve the robustness of the algorithm, especially with strongly varying terms (gravity affects very differently a denser zone with respect to a less dense one).
In other words, it is a way to improve stability when enforcing the force balance.

Best,

[mdjames]

Oh, okay. I think the first post makes sense to me, if i understand the new variables. The second one definitely does.

Thanks again...you're right, I would like to know the correct way to do it rather than a way.

Next up is getting my forcing term time-dependent..."here be dragons"

I'll check back if I can't manage to work something out.

[alberto]

Originally Posted by mdjames

Oh, okay. I think the first post makes sense to me, if i understand the new variables. The second one definitely does.

Ask if something is not clear :-)

Thanks again...you're right, I would like to know the correct way to do it rather than a way.

Let's say nothing is written in stone. If you work with incompressible single-phase codes, you probably won't see any difference with the two approaches. Things are very different in compressible or multiphase flows.

Next up is getting my forcing term time-dependent..."here be dragons"

I'll check back if I can't manage to work something out.

The principle is the same, but of course ask if you have problems.

Best,

[mdjames]

I think I'm hot on the trail, I just can't find to seem to find a way to call out the runtime as a dimensioned scalar. I tried creating a new time variable with:

dimensionedScalar timemdj
(
timemdj [0 0 1 0 0 0 0 ] runTime.Value()
);

to include in my equation for the body force but to no avail...

This is getting fun!

[alberto]

It is

dimensionedScalar timemdj
(
"timemdj",
dimensionSet(0, 0, 1, 0, 0),
<yourExpression> // Replace with what you need
);

or, simply

dimensionedScalar timemdj = <yourExpression>;

Best,

[mdjames]

I'm confused...

my forcing is described as

dimensionedVector waveP=amp*cos(-freq*time_); //determining the new forcing

and this yields a forcing in the x-direction as I wanted but I realized:

1. I never included any information to make this an x forcing
2. I'm adding a vector to a matrix (in the UEqn, below)???

tmp<fvVectorMatrix> UEqn
(
fvm::ddt(U)
+ fvm::div(phi, U)
+ turbulence->divDevReff(U)
- waveP
);

Is there some assumption on OF's part or something?

[alberto]

How do you define Amp? It has to be a vector for the code to compile...

Best,

[mdjames]

Okay, I see now...minor freak-out.

I don't quite understand how the vector waveP is added to higher dimension arrays like U...

[alberto]

Originally Posted by mdjames

Okay, I see now...minor freak-out.

I don't quite understand how the vector waveP is added to higher dimension arrays like U...

waveP is a vector of three components, since amp has to be a vector of three components, multiplied by the time-dependent part.

Each component of waveP is added to the corresponding scalar equation for the component of U. So, say you have

g= (gx, gy, gz)

and you do

UEqn += g;

You add gx to the equation for Ux, gy to the equation for Uy and so on. OpenFOAM syntax just makes it easier for you to do this in only one step.

Best,

[mdjames]

Okay, to check my understanding of the grand scheme of things:

Each node has a vector equation associated with it. So, in effect we have the matrix of equations (UEqn) added to a matrix of the same size, with all elements being the vector waveP

correct?

If so, I suppose creating a spatially variable forcing would just involve directly specifying a forcing matrix.

Am I on the right track here?

Thanks a ton!

[alberto]

Yes. Put into the code, you simply define a volVectorField, initialize it with the values of your force term, and add it to UEqn.

P.S. You amp is (|amp|, 0, 0) right?

[mdjames]

Beautiful! Seems like this may be a fairly painless solution. Much simpler than some other schemes I had rolling around in my head.

Yep, just without the commas It's read from a wavePressDict like viscosity is.

[Thomas Baumann]

Hi,

I'm trying to implement new anisotropic turbulence models in OpenFOAM.
It's like a normal isotropic turbulence-model, with an explicit nonlinear part too.

So the Reynolds-fluxes are defined as: R = linearStress + nonlinearStress (explicit)

so divDevReff=
- fvm::laplacian(nuEff(), U)
- fvc::div(nuEff()*dev(fvc::grad(U)().T()))
+ fvc::div(nonlinearStress_)

Would it make here sense to put the nonlinear-part in the flux-part as descirbed above? Would the results be better than putting all of them in a fvVectorMatrix? Or is this only for the convergence, the robustness of the system better?

There are some nonlinear models in OF-1.7 like the nonlinearshih model, but here the nonlinearStress_-term is in the fvVectorMatrix, too. Using them to solve a fully developed channel flow I have velocities out of the wall. Could it be the reason, because the nonlinearStress-part is described in the fvVectorMatrix?

Thanks and regards,
Thomas

[be_inspired]

Hi Alberto,

In case the force is managed in the first way ( force term in the flux), my question is if Ueqn must be modified or not.

I am working over the rotorDiskSource (OF2.3.0) and some pressure wiggles happen when discreted body forces are included directly over the Ueqn.
I think that your approach could improve this situation.
http://www.openfoam.org/mantisbt/view.php?id=1313

A*U = H + F

I think it should be:
A*U = H - grad[p] + F
U=H/A - 1/A*grad[p] + F/A

Excuse me for using the same topic but there are other topics where others have reported this bug and no answer have been given.

http://www.cfd-online.com/Forums/ope...rce-model.html
http://www.cfd-online.com/Forums/ope...l-wiggles.html

I am trying to implement the body force in this other way and I am not able to do it correctly. Could you please help me?

[manoj jaiswal]

Hi i am a new Foamer and i am simulating flow over a cylinder(in pimplefoam) . i want to add a body force in the flow , so what are the steps needed for the same . I have gone through a tutorial of how to add temperature in icoFoam (http://openfoamwiki.net/index.php/Ho...ure_to_icoFoam) but still not clear about the addition of body force.

Thanks

Originally Posted by alberto

Including force terms in the flux is a way to improve the robustness of the algorithm in some cases (for example when the force is related to the density, like for the gravity).

If you write your momentum equation as

ddt(U) + div(UU) = div(tau) - grad(p) + F

where F is your forcing term, and then you write in semi-discrete form you get

A*U = H + F

where I put everything I do not want to carry around in H. Now, the predicted velocity is

U = H/A + F/A

and the flux, obtained interpolating this on the cell faces (indicated by _f) is

phi = (H/A)_f . S + (F/A)_f . S

In your case F is constant, so

phi = (H/A)_f . S + F(1/A)_f . S

In OpenFOAM notation, this means that you should add

surfaceScalarField rUAf = fvc::interpolate(rUAf)

surfaceScalarField phiForcing = rUAf*(F & mesh.Sf())

to your flux, after adjustPhi. At this point however you have to modify the correction step for the velocity too (the flux is fine), consistently with the new velocity predictor:

U += fvc::reconstruct(phiForcing) - rUA*fvc::grad(p);

On a side note, you probably won't need this in your case, but it is better to know

Hopefully I did not introduce errors while typing

Best,