[asaha]

Hello Armin,

I hope the link will be of any help.
http://openfoamwiki.net/images/4/4a/Simple2DFlow.tar.gz

[armin_h]

Dear Saha

Mnay Rhanks about the comments.
Appreciate about that.

Reagrds
Armin

[ehsanyasari]

Hello all
I am going to run a channel395 with channel oodles solver, when i choose mixedsmagorinskey model and laplace filter, after some iteration, courant number diverged. but when i used simple filter it works,and also when i changed the mesh into uniform mesh, it works.


I would appreciate any comments.
Best Regards
Ehsan

[grtabor]

Dear All,

I am running (or attempting to run at any rate) a surface-mounted cube case using LES. I have it running successfully on a fairly coarse mesh (72000 cells) with a timestep which gives a maximum Courant number of around 0.1. I attempted to refine the mesh somewhat, taking it to 140,000 cells, then mapping the existing solution across and restarting with a timestep 1/5th of the original. Trouble is, after about 20 timesteps (with max courant number of around 0.06) the solution blows up, and I cannot make sense of what the problem is. Initially the k field seemed to be going unbounded, followed by the Courant number exploding; however swapping SGS models (from oneEqEddy to Smagorinsky or to dynOneEqEddy) makes no difference. I have also tried drastically smaller timesteps, non-orthogonal correctors, and some changes to the differencing schemes used; without success. Any other suggestions as to what might be amiss here would be gratefully accepted.

On a broader note; is there any source of information about what is available in the LES modelling area in OpenFOAM at the moment? The Users Guide is rather terse on this issue, and the usual reference papers from the 1990's are possibly out of date. I'd be particularly interested in near-wall modelling.

Gavin

[andrewburns]

What does checkMesh say about your new finer mesh? I haven't got any experience with LES but generally when my solution blows up with unbounding k and epsilon it's due to a poor mesh. You could try switching from linear corrected to linear limited 0.5 or so for k and epsilon too.

[braennstroem]

Hi,

I would like to take a look at my mesh resolution and check, if it is fine enough.

There are actually a couple of different methods like the checking the energy spectra .

Another way is the to check the turbulence resolution, i.e. the ratio of modelled to resolved turbulent kinetic energy. Does anyone know, how I can postprocess this using openfoam?

Fabian

[grtabor]

Sounds fairly straightforward. You have the U and k fields; write a code to read these in and calculate 0.5*sqr(U)/k.

Gavin

[maka]

@(August 10, 2005 - 04:50 pm)
\quote(nuSgs doesn't matter because it is calculated from k and delta.)

I have a related question. I've been using channelOodles for a while now. But suddenly, I asked myself why do we need to specify an initial and boundary condition for nuSgs if it is calculated from resolved (like Sij) and modeled variables like k) variables which have an initial and boundary condition of their own.

I know that in case of wall functions for example specifying a boundary condition on nuSgs is useful but will such b.c. be consistent with the suface value that can be calculated from the resolved and modeld fields. I'm confused :-(. Would you please give me some feedback on this.

Best regards,
Maka.

[braennstroem]

Hi Gavin,

thanks, fairly straightforward; I would do it as a postprocessing step, but using the Smagorinsky I get no written k for the internal field. I actually tried to add:

k_
(
IOobject
(
"k",
runTime_.timeName(),
mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh_
),

and added:

volScalarField k_;

to the Smagorinsky class...

but this is not enough. Do you have an idea?

Fabian

[grtabor]

Hi Fabian,

You are right - Smag does not use k. I'm tempted to suggest you switch to the 1-equation model, at least for a bit to check the mesh. Alternatively, the basic equation for the viscosity from the Smag model is

nu = C_s^2 \delta^2 |S|

and you have nu, so I guess you could fiddle with this to get an estimate for k. Was this what you tried? I wouldn't bother adding it to the smagorinsky class - overkill, and it breaks the elegant relationship between the mathematical structure of the various models and the class inheretance (sorry; I designed that bit of the code!)

Gavin

[braennstroem]

Hi Gavin,

Thanks!
Actually the implementation of Smagorinsky uses:

//- Return SGS kinetic energy
// calculated from the given velocity gradient
tmp<volscalarfield> k(const tmp<voltensorfield>& gradU) const
{
return (2.0*ck_/ce_)*sqr(delta())*magSqr(dev(symm(gradU)));
}

//- Return SGS kinetic energy
tmp<volscalarfield> k() const
{
return k(fvc::grad(U()));
}

And I thought I could write these values to a file using the earlier mentioned lines for k!?

Fabian

[ville]

Hi,
what do you have in the turbulenceProperties dictionary as the LESmodel?
-Ville

[braennstroem]

Hi Ville,
I know 'turbulenceProperties' just in combination with the actual calculation!? The above problem occurs when compiling a small post-tool, which needs the filter width.

Fabian

[braennstroem]

I actually solved the above problem using these lines:


#include "fvCFD.H"
#include "incompressible/LESmodel/LESmodel.H"
#include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
#include "incompressible/transportModel/transportModel.H"
dimensionedScalar smallU ("smallU", dimensionSet(0,1,-1,0,0), SMALL);

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



int main(int argc, char *argv[])
{

# include "addTimeOptions.H"
# include "setRootCase.H"

# include "createTime.H"

// Get times list
instantList Times = runTime.times();

// set startTime and endTime depending on -time and -latestTime options
# include "checkTimeOptions.H"

runTime.setTime(Times[startTime], startTime);

# include "createMesh.H"
# include "turbRes.H"
Info<< "delta = " <<
sgsModel->delta() <<endl;

return(0);
}

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



and the 'turbRes.H' file:



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
);

# include "createPhi.H"


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


singlePhaseTransportModel laminarTransport(U, phi);

autoPtr<lesmodel> sgsModel
(
LESmodel::New(U, phi, laminarTransport)
);


I get problems, when running the code. It tries to open a 'k' file, which is not available, but where in these lines does this happen!?

Fabian

[ngj]

Hi all

Before going around submitting it as a bug-report, I'd like to discuss if I am off track first.

In nuSgsWallFunctionFvPatchScalarField.C the friction velocity is calculated as:

scalar utau = sqrt((nuSgsw[facei] + nuw[facei])*magFaceGradU[facei]);

with magFaceGradU being U.snGrad().

In anything but steady state and exceptionally simple geometries, there will be a normal velocity component at the wall, thus the gradient in this normal to the wall is incorporated in U.snGrad(), thus

mag(U.snGrad()) != mag(Uf).

Shouldn't U.snGrad() be projected onto the wall to avoid the normal stress in the bed shear stress?

I am looking forward to hear your point of view.

Best regards,

Niels

[eugene]

This issue is known. The rationalization for keeping it as it is is as follows:

1. The wall normal gradient of velocity will tend to be small compared to the tangential gradient.
2. If there is a significant wall normal gradient then the wall functions are no longer valid in any case and the relative contribution of the error made by including the magnitude of the wall normal gradient is diminished due to the error incurred by assuming equilibrium.

All that said, subtracting the normal component before calculating the gradient would be more correct.

[ngj]

Thanks for the explanation.

Would you have any reference on the implemented wall function?

Best regards,

Niels

[maka]

[vanDriestDelta]
What does wallPointYPlus::yPlusCutOff = 500; do? any reason behind that? Thanks.

Best regards,
Maka.

[ngj]

Hi Maka

In the vanDriest, a damping function given as exp(-yPlus) is used (cannot recall if that is the exact expression), anyway, as yPlus exceeds some value, exp(-yPlus) is effectively 0.

The cutOff is used in the this part of the code:

/OpenFOAM/OpenFOAM-1.4.1/src/finiteVolume/fvMesh/wallDist/wallPointYPlus

as this calculates the distance to the wall in yPlus-coordinates, it also needs to know the bed shear stress in the nearest wall face, thus as I see it, it is a upper limit for this interconnection between distance and wall.

Best regards,

Niels

[maka]

The expression is:

delta_ = min
(
static_cast<const>(geometricDelta_()),
(kappa_/Cdelta_)*((scalar(1) + SMALL) - exp(-y/ystar/Aplus_))*y
);

at high y/ystar the exp term is zero and the expression reduces to:
min(geometricDelta_(), (kappa_/Cdelta_)*y)

I tried to understand what the cutOff do to y from
src/finiteVolume/fvMesh/wallDist/wallPointYPlus

but it was not clear. What I understood is that after the cutOff value it does not update y?

// only propagate if interesting (i.e. y+ < 100)
scalar yPlus = Foam::sqrt(dist2)/w2.data();


if (yPlus < yPlusCutOff)
{
// update with new values
distSqr() = dist2;
origin() = w2.origin();
data() = w2.data();

return true;
}

But what is the value of y when it is not updated? Thanks for your help.