Hi all,
I am a beginner in OpenFoam. I am trying to simulate waves and as a convergence and stability test I try to run the simple case of a standing wave in a 1x1 tank.
I tried to increase the Mesh accuracy and I measure the convergence two ways :
with the Kinematic energy : I measure it at each time step by modifying the solver .C file (ihFoam in my case but it should give the same result with interFoam) by adding the lines :
Code:
OFstream energyTotIntegral("energyTotIntegral.dat");
while (runTime.run())
{
scalar QQ = fvc::domainIntegrate(rho*mag(U)*mag(U)*.5*alpha1).value();
energyCinIntegral << runTime.timeName() << " "<< QQ << endl;
and a L2 convergence test on the velocity profile :
Code:
double ux (double a, double om, double k, double h, double time, double x, double z)
{
double res =0;
res = a*om*Foam::cosh(k*(h+z))/Foam::sinh(k*h);
res = res*Foam::sin(k*x)*Foam::sin(om*time);
// Info<< x << " " << z << " "<< time <<" " << res << endl;
return res;
}
to calculate the standing waves Airy (Stokes I) solution for small perturbations and I measure the error like that :
Code:
// Calculate the error in velocity
//First Construct the volscalar field of Ux-uAnalytic
volScalarField Ux = U.component(1);
volScalarField Uz = U.component(2);
volScalarField xcoord = mesh.C().component(1);
volScalarField zcoord = mesh.C().component(2);
// Info<< "Here" << endl;
forAll (mesh.C(), celli)
{
// Info<< "Here" << endl;
Ux[celli]=alpha1[celli]*(Ux[celli]-ux(a,om,k,h,runTime.value(),xcoord[celli],zcoord[celli]));
Ux[celli]=Ux[celli]*Ux[celli];
Uz[celli]=alpha1[celli]*(Uz[celli]-uz(a,om,k,h,runTime.value(),xcoord[celli],zcoord[celli]));
Uz[celli]=Uz[celli]*Uz[celli];
}
scalar er1 = fvc::domainIntegrate(Ux).value();
errorUx << runTime.timeName() << " "<< er1 << endl;
scalar er2 = fvc::domainIntegrate(Uz).value();
errorUz << runTime.timeName() << " "<< er2 << endl;
My mesh is the following :
Code:
convertToMeters 1;
vertices
(
(0.500 -0.5 -1 )
(0.505 -0.5 -1 )
(0.505 0.5 -1 )
(0.500 0.5 -1 )
(0.500 -0.5 0.2 )
(0.505 -0.5 0.2 )
(0.505 0.5 0.2 )
(0.500 0.5 0.2 )
);
blocks
(
hex (0 1 2 3 4 5 6 7) (1 100 120) simpleGrading (1 1 1)
);
The BC are :
Code:
patches
(
patch atmosphere
(
(5 6 7 4)
)
empty frontAndBack
(
//Back
(0 4 7 3)
//Front
(1 2 6 5)
)
wall defaultFaces
(
)
);
and I set up
for alpha1 :
Code:
dimensions [0 0 0 0 0 0 0];
internalField uniform 0;
boundaryField
{
atmosphere
{
type inletOutlet;
inletValue uniform 0;
value uniform 0;
}
defaultFaces
{
type zeroGradient;
}
frontAndBack
{
type empty;
}
}
for p_rgh
Code:
dimensions [1 -1 -2 0 0 0 0];
internalField uniform 0;
boundaryField
{
atmosphere
{
//type fixedValue;
//value uniform 0;
type totalPressure;
U U;
phi phi;
rho rho;
psi none;
gamma 1;
p0 uniform 0;
value uniform 0;
}
defaultFaces
{
type slip;
}
frontAndBack
{
type empty;
}
and for U
Code:
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
atmosphere
{
type pressureInletOutletVelocity;
value uniform (0 0 0);
}
defaultFaces
{
type slip;
}
frontAndBack
{
type empty;
}
}
Finally the schemes I use are the following :
Code:
ddtSchemes
{
default CrankNicolson/Euler or backward;
}
gradSchemes
{
default cellMDLimited Gauss linear 1;
}
divSchemes
{
div((rho*phi|interpolate(porosity)),U) Gauss limitedLinearV 1;
div(phi,alpha) Gauss vanLeer01;
div(phirb,alpha) Gauss interfaceCompression;
}
laplacianSchemes
{
default Gauss linear corrected;
}
interpolationSchemes
{
default linear;
}
snGradSchemes
{
default corrected;
}
fluxRequired
{
default no;
p_rgh;
pcorr;
alpha1;
}
Ok now Initial Conditions : I set up the field using the command :
Code:
echo Setting the fields...
funkySetFields -time 0 -field alpha1 -expression '0.5-tanh((pos().z-(0.+0.02*cos(2*pi*(pos().y))))/(3/4*1/N))*0.5' -time 0 > funkySetField.log
the tanh should spread the interface on 3 cells and avoid "loss" of mass at the beginning of the simulation.
All the rest of the fields are set to zero.
So I have :
h=1m
A=0.02 m << h
L=1
Dispersion relation gives om = 7.85096 rad/s so f = 1.24952 Hz
I run the simulations for 7 seconds
and I tried the cases
dt = {0.01,0.005, 0.001}
and L/Dy = L/Dz = {50,100,200,500}=N
I use the scheme Euler first and I get rather good results. It is stable for all the cases but dt=0.01 and N=500.
Unfortunately the case backward does not work...
I played with the case CrankNicolson and it is unstable if I put a too big Psi.
Can somebody help me out? Where am I wrong? Wrong BCs? Wrong Schemes?
If you are interested I can send you the case by email for let's say N=100 and dt =0.005 so it is faster for you to try the case.
Thanks
Remi
