CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (http://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Approximate Deconvolution model for wall bounded flows channeloodles (http://www.cfd-online.com/Forums/openfoam-solving/60006-approximate-deconvolution-model-wall-bounded-flows-channeloodles.html)

turnow September 19, 2006 04:00

Hi all, I am currently tryi
 
Hi all,

I am currently trying to implement a Model called "Approximate deconvolution model (ADM)" for wall bounded flows. The ADM model for large eddy simulations has been proposed by S. Stolz and N.Adams. The main ingredient is an approximation of the nonfiltered Field by truncated series expansion of the inverse Filter operator.

Therefore the momentum equation in the channeloodles solver becomes explicit.
After running the code for 10 - 15 timesteps (deltaT = 0.5e-03) the simulation becomes instable and i get pressure peaks in the corners of the discretized channel.

Here is modified channelOodles solver
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //


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

for(runTime++; !runTime.end(); runTime++)
{
Info<< "Time = " << runTime.timeName() << nl << endl;


# include "CourantNo.H"

# include "readPISOControls.H"

//sgsModel->correct();

autoPtr<lesfilter> gridFilterPtr_
(
LESfilter::New(U.mesh(),sgsModel->LESmodelProperties())
);

LESfilter& gridFilter_(gridFilterPtr_());
dimensionedScalar Xsi
(
"Xsi",
dimensionSet(0, 0, -1, 0, 0),
1.0
);

//the deconvolved field

Ustar = 3.0*(U) - 3.0*(gridFilter_(U)) + gridFilter_(gridFilter_(U));

fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvc::div(gridFilter_(Ustar * Ustar)(), "div(phi,U)")
- fvc::div(nu*dev(symm(fvc::grad(U))), "div(diss,U)")
- Xsi*(gridFilter_(Ustar) - U)
==
flowDirection*gradP
);

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

# include "continuityErrs.H"


// --- PISO

volScalarField rUA = 1.0/UEqn.A();

for (int corr=0; corr<nCorr; corr++)
{
U = rUA*UEqn.H();
phi = (fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, U, phi);

for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(rUA,p) == fvc::div(phi)
);

pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();

if (nonOrth == nNonOrthCorr)
{
phi -= pEqn.flux();
}
}

# include "continuityErrs.H"

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

U.correctBoundaryConditions();
}


Please don't look at the programstyle, I am still working on it.

Thanks for your help.
J. Turnow

david_h September 19, 2006 10:06

Have you tried adding some imp
 
Have you tried adding some implicitness to the momentum equation ? i.e. "fvm" instead of "fvc" for the viscous term and replacing the relaxation term:
- Xsi*(gridFilter_(Ustar) - U)
with
-Xsi*gridFilter_(Ustar) + fvm::Sp(Xsi,U)

turnow September 19, 2006 10:55

Hi for the visous term "gr
 
Hi

for the visous term
"grad(U)" is not possible for fvm:: ?
Maybe you have an idea to avoid using fvm::grad(U)?

Furthermore I corrected the convective term in the momentum equation to

+ fvc::div(fvc::interpolate(gridFilter_(Ustar * Ustar)) & mesh.Sf())

I also tried a new explicit pressure correction,
but I got still pressure peaks at the corners of the channel.


volVectorField flux = fvc::div(fvc::interpolate(gridFilter_(Ustar * Ustar)) & mesh.Sf());

for (int corr=0; corr<nCorr; corr++)
{
fvScalarMatrix pEqn
(
fvm::laplacian(p) == -fvc::div(flux + Xsi*(U - gridFilter_(Ustar)))
);

pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();
}

Thanks for all your help.
J. Turnow

eugene September 19, 2006 12:09

Are you using a plane channel
 
Are you using a plane channel test case? If so, try changing the side boundaries to walls.

What happens if you make your timestep 10 times smaller?

I can't see anything obviously wrong, although I must admit I don't know the method.

david_h September 19, 2006 13:41

You are quite right about fvm:
 
You are quite right about fvm::grad.

I meant something along lines of the following:
fvm::laplacian(nu,U) + fvc::div(nu*dev(fvc::grad(U)) )
which is similar to the code for the divR() function in "turbulenceModels/incompressible/laminar/laminar.C"

Dave

turnow September 20, 2006 03:07

thanks eugene and dave, Yes
 
thanks eugene and dave,

Yes I am using the plane channel test case.
I changed the side boundary's to walls.
But nothing happend. Again my pressure values blow up to -300 in the first 2 timesteps.
By reducing the time step to 0.5e-04 it is the same game.

Also using the implicitness for the viscous term
"- fvm::laplacian(nu, U) - fvc::div(nu*dev(fvc::grad(U)().T()))"
the simulations is getting instable.

I am wondering also about the term

pEqn.setReference(pRefCell, pRefValue);

When using the term my pressure is blowing up to 300, but by neglecting the reference term the pressure rises up to -3.47146e+13.
Another point is the velocity field. Here it is not effected by the high pressure values.
In addition I compared the convective terms from the ADM solver and the original channelOodles solver. The values are in the same range.

If you have access the Journal "Physics of Fluids" and you would like to read more about the ADM model, there's a good articel:
"An approximate deconvolution model for large-eddy simulation with application to incompressibe wall bounded flows" ( Volume 13, Number 4, April 2001// page 997-1015)

eugene September 20, 2006 05:49

Hmm, I seem to recall there wa
 
Hmm, I seem to recall there was an issue with the reference cell being located adjacent to a cyclic boundary. Try moving the reference cell to somewhere in the middle of the domain.

On the other hand, the fact that it happens at all the corners seems to indicate that the origin of the issue is elsewhere.

If the above doesn't work, you will have to start stripping stuff out until you are left with the minimum system that reproduces the problem. The first step would be to rewrite the standard channelOodles explicitly to see if you get the same issue.


All times are GMT -4. The time now is 18:08.