[spaceprop]

Setup:
~15cm diameter x ~30cm long cylindrical tank, water, 0 gravity, surface tension, 60deg contact angle.
~800k cell mesh, hex dominant, 6 prism layers on walls, good transition from prism layer to hex layer, generated with SHM, checkMesh approved.
Sinusoidal excitation: 1 axis, 10cm amplitude, 0.25Hz, 1000 position points per second, 10s runtime.
Standard force calculation utility.
Settings based on sloshingtank tutorials, see below.
Hierarchical decomposition
Parallel interDyMFoam run on 16 cores (1 node).

Description of problem:
The case seems to run fine. All parameters converge (initial residuals) to at least 1e-4 within 10-20 pimple iterations. The fluid interface at all time steps looks reasonable, aka, nothing seems to be blowing up. HOWEVER, the axial force waveform shows very high amplitude and high frequency (on order of 1/time step) oscillations, most of which comes from the pressure component. Attached are pictures of the total force (p+visc) along the excitation direction. One is unfiltered and one is filtered with a 20Hz cutoff, 4th order, Cheby II low pass filter. I think the filtered one may be somewhat representative of what is actually happening. The weird thing is, I don't recall seeing these oscillations in the sloshingtank tutorial cases (with 1g gravity).

1. Has anyone seen this oscillation phenomenon before?
2. Any ideas why it happens? My guess is that it has something to do with fluid surface interface tracking, but I honestly don't know.

Other unrelated questions:
3. Anyone know if you can have a pimple residual control or relaxation factor for alpha.water? (commented out in my fvSolution)
4. For this case, is p_rgh or p the correct field for the relaxation factor?

fvSolution:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  3.0.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    alpha.water
    {
        nAlphaCorr      3;
        nAlphaSubCycles 1; //must be 1 for crank-nicolson
        cAlpha          1;

//might make more accurate:	
	MULESCorr       yes;
        nLimiterIter    10;
  	alphaApplyPrevCorr  no;

        solver          smoothSolver;
        smoother        symGaussSeidel;
        tolerance       1e-9;
        relTol          0;
    }
    alpha.waterFinal
    {
	cAlpha          1;	
	solver          smoothSolver;
        smoother        symGaussSeidel;
	tolerance       2e-10;
        relTol          0;
    }

    pcorr
    {
        solver          PCG;
        preconditioner
        {
            preconditioner  GAMG;
            tolerance       1e-05;
            relTol          0;
            smoother        DICGaussSeidel;
            nPreSweeps      0;
            nPostSweeps     2;
            nFinestSweeps   2;
            cacheAgglomeration false;
            nCellsInCoarsestLevel 10;
            agglomerator    faceAreaPair;
            mergeLevels     1;
        }

        tolerance       1e-05;
        relTol          0;
        maxIter         100;
    }

    p_rgh
    {
        solver          GAMG;
        tolerance       1e-08;
        relTol          0.001;
        smoother        DICGaussSeidel;
        nPreSweeps      0;
        nPostSweeps     2;
        nFinestSweeps   2;
        cacheAgglomeration true;
        nCellsInCoarsestLevel 128;
        agglomerator    faceAreaPair;
        mergeLevels     1;
        maxIter		100;
    }

    p_rghFinal
    {
        solver          PCG;
        preconditioner
        {
            preconditioner  GAMG;
            tolerance       2e-09;
            relTol          0;
            nVcycles        2;
            smoother        DICGaussSeidel;
            nPreSweeps      2;
            nPostSweeps     2;
            nFinestSweeps   2;
            cacheAgglomeration true;
            nCellsInCoarsestLevel 128;
            agglomerator    faceAreaPair;
            mergeLevels     1;
        }

        tolerance       2e-09;
        relTol          0;
        maxIter         20;
    }

    U
    {
        solver          smoothSolver;
        smoother        GaussSeidel;
        tolerance       1e-07;
        relTol          0;
        nSweeps         1;
    }

    UFinal
    {
	$U;
        tolerance       1e-08;
        relTol          0;
    }
}

PIMPLE
{
    momentumPredictor yes; //solving this might be more accurate. 
    nOuterCorrectors 20;
    nCorrectors     2;
    nNonOrthogonalCorrectors 1;
    correctPhi      no;

    pRefPoint       (0 0 0.15);
    pRefValue       1e5;

    residualControl
    {

	p_rgh
	{
	    tolerance 1e-4;
	    relTol 0;
	}
	U
	{
	    tolerance 1e-4;
	    relTol 0;
	}
	//alpha.water
	//{
	//    tolerance 1e-3;
	//    relTol 0;
	//}

    }

}

relaxationFactors
{
    
    fields
    {
	p_rgh		0.3;
	p		0.3;
	p_rghFinal	1;
	pFinal		1;
	//alpha.water	0.7;
    }
    equations
    {
        UFinal	1;
	U	0.6;
    }
}


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

fvScheme:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  3.0.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         CrankNicolson 0.9; //was Euler (1st order)
}

gradSchemes
{
    default         Gauss linear;
}

divSchemes
{
    div(rhoPhi,U)  Gauss linearUpwind grad(U); //was vanLeerV. lots of options
    div(phi,alpha)  Gauss vanLeer; //can only be vanLeer
    div(phirb,alpha) Gauss vanLeer; //linear works. was vanLeer
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear; //can only be linear
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

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

controlDict:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  3.0.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      controlDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

application     interDyMFoam;

startFrom       startTime;

startTime       0;

stopAt          endTime;

endTime         10;

deltaT          0.001;

writeControl    adjustableRunTime;

writeInterval   0.1;

purgeWrite      0;

writeFormat     ascii;

writePrecision  7;

writeCompression compressed;

timeFormat      general;

timePrecision   6;

runTimeModifiable yes;

adjustTimeStep  yes;

maxCo           1.5;
maxAlphaCo      1.5;
maxDeltaT       0.01;

functions
{
    #include "forcesIncompressible";
    #include "surfaces";
    //#include "residuals";
}


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

forcesIncompressible:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  3.0.1                                 |  
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      forcesIncompressible;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

forces
{
    rhoInf      997.561;    // Fluid density
    patches     ( walls);

    CofR        (0 0 0);
    pitchAxis   (0 1 0);

    #include "forces.cfg"
}

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

forces.cfg:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  3.0.1                                 |  
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    object      forces.cfg;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

type            forces;
functionObjectLibs ( "libforces.so" );

enabled         true;
outputControl   timeStep;
outputInterval  1;

pName           p;
UName           U;
rhoName         rhoInf; // Incompressible solver
log             on;

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

[spaceprop]

Running 2.4M and 6M cell cases now. If anything the amplitude of oscillation increases with increasing mesh size.

The number of pimple iterations required for 1e-4 convergence also increased to 15-25 for the larger mesh cases. The Uz (motion is in z direction) or p equation tends to have slower convergence.

[spaceprop]

I've made some progress on this.

Short:
The default precision used in writing the 6DoF.dat file (from gen6DoF.c) is not high enough: 6 digits is stream operator default. While this may produce a smooth looking position curve, any numerical imperfections are amplified by double differentiation to acceleration to form high frequency, high amplitude oscillations. Since F=m*a, the force calculations are picking up this numerical "noise". However, I don't think this accounts for all of the oscillations.

Details:
Test case:
gen6DoF.c: modified the translation equation to be (1-cos(f*t))*A so that position and velocity start at 0. Z-axis motion only: 10cm amplitude, 0.25Hz, 1000 points per second, 10s.

Results:
If I numerically differentiate (central diff) the position data to get velocity, and then numerically differentiate velocity to get acceleration, I see the "Default" graph attached. To get the "MATLAB" plot, I mimicked the routine in gen6DoF.c in MATLAB, wrote the output file with 12 decimal point precision, then read it back in and repeated the numerical differentiation. The result is a smooth acceleration curve. Note: 8 decimal point precision was not enough. High frequency oscillations were still visible in the acceleration plot, though their amplitude was significantly smaller.

Recommendations:
If you are using gen6DoF for anything requiring precision motion and/or force measurements, I strongly suggest increasing the precision (put setprecision(x) in stream, where x is at least 12). I will be using MATLAB because I have trajectories outside the scope of gen6DoF to create. If anyone is interested in the MATLAB gen6DoF code, here ya go:

ts=0.001;
endtime=10;
ftrans=[0 0 0.25]; %[x y z] Hz
fomega=[0 0 0]; %[x y z] rad/s
Pi=3.14159265359;
Atrans=[0 0 0.1]; %[x y z] m
Aomega=[0 0 0]; %[x y z] deg

t=0:ts:endtime;
lent=length(t);
POS=(1-cos(2*Pi*t'*ftrans)).*repmat(Atrans,lent,1);
ANG=sin(t'*fomega).*repmat(Aomega,lent,1);
% xpos=POS(:,1);
% ypos=POS(:,2);
% zpos=POS(:,3);
% vel=gradient(zpos,ts);
% acc=gradient(vel,ts);

fileID = fopen('6DoF2.dat','w');
fprintf(fileID,'%d\n',lent);
fprintf(fileID,'%s\n','(');
for i=1:1:lent
    fprintf(fileID,'%s%.12g%s%.12g %.12g %.12g%s%.12g %.12g %.12g%s\n','(',...
        t(i),' ((',POS(i,1),POS(i,2),POS(i,3),') (',ANG(i,1),ANG(i,2),ANG(i,3),')))');
end
fprintf(fileID,'%s\n',')');
fclose(fileID);

Simulation progress:
My simulations are showing significantly less high frequency noise in the force waveforms, though some is still there (about 2 orders of magnitude lower amplitude). I will post plots when the simulations have progressed farther. My guesses as to what this could be due to:
1. Perhaps 12 decimal precision is not enough, though double precision has only 15-17 significant digits of decimal precision.
2. Numerical noise from other calculations.
3. Interface tracking scheme is slightly oscillatory. There are more advanced (and computationally expensive) schemes available.

I did another test where I ran a transient no-motion simulation (interFoam) with the 800k cell mesh until the surface seemed mostly steady. Some motion of the fluid interface, but never more than 1 cell, and no waves/oscillations. See attached plot of the forces. I think those oscillations are probably due to the interface tracking scheme.

Thoughts anybody?

[spaceprop]

Another related question: What is the precision of the values being read from the 6dof.dat file?

I spent an hour or so digging through the source code. I found the tabulated6dofmotion source, and the interpolatesplinexy source, and that lead me to the scalarfield source. And that lists either single precision floats or double precision floats field types. But what controls which type, SP or DP, the vector components are? I know that if you declare a variable as double, then use an ifstream, it will read in as double, but what in openfoam is doing the declaration? I compiled with DP flag on: does that mean that all variables in openfoam, including anything being read in from tables, are read in as doubles?

[Di Mu]

Originally Posted by spaceprop

I've made some progress on this.

Short:
The default precision used in writing the 6DoF.dat file (from gen6DoF.c) is not high enough: 6 digits is stream operator default. While this may produce a smooth looking position curve, any numerical imperfections are amplified by double differentiation to acceleration to form high frequency, high amplitude oscillations. Since F=m*a, the force calculations are picking up this numerical "noise". However, I don't think this accounts for all of the oscillations.

Details:
Test case:
gen6DoF.c: modified the translation equation to be (1-cos(f*t))*A so that position and velocity start at 0. Z-axis motion only: 10cm amplitude, 0.25Hz, 1000 points per second, 10s.

Results:
If I numerically differentiate (central diff) the position data to get velocity, and then numerically differentiate velocity to get acceleration, I see the "Default" graph attached. To get the "MATLAB" plot, I mimicked the routine in gen6DoF.c in MATLAB, wrote the output file with 12 decimal point precision, then read it back in and repeated the numerical differentiation. The result is a smooth acceleration curve. Note: 8 decimal point precision was not enough. High frequency oscillations were still visible in the acceleration plot, though their amplitude was significantly smaller.

Recommendations:
If you are using gen6DoF for anything requiring precision motion and/or force measurements, I strongly suggest increasing the precision (put setprecision(x) in stream, where x is at least 12). I will be using MATLAB because I have trajectories outside the scope of gen6DoF to create. If anyone is interested in the MATLAB gen6DoF code, here ya go:

ts=0.001;
endtime=10;
ftrans=[0 0 0.25]; %[x y z] Hz
fomega=[0 0 0]; %[x y z] rad/s
Pi=3.14159265359;
Atrans=[0 0 0.1]; %[x y z] m
Aomega=[0 0 0]; %[x y z] deg

t=0:ts:endtime;
lent=length(t);
POS=(1-cos(2*Pi*t'*ftrans)).*repmat(Atrans,lent,1);
ANG=sin(t'*fomega).*repmat(Aomega,lent,1);
% xpos=POS(:,1);
% ypos=POS(:,2);
% zpos=POS(:,3);
% vel=gradient(zpos,ts);
% acc=gradient(vel,ts);

fileID = fopen('6DoF2.dat','w');
fprintf(fileID,'%d\n',lent);
fprintf(fileID,'%s\n','(');
for i=1:1:lent
    fprintf(fileID,'%s%.12g%s%.12g %.12g %.12g%s%.12g %.12g %.12g%s\n','(',...
        t(i),' ((',POS(i,1),POS(i,2),POS(i,3),') (',ANG(i,1),ANG(i,2),ANG(i,3),')))');
end
fprintf(fileID,'%s\n',')');
fclose(fileID);

Simulation progress:
My simulations are showing significantly less high frequency noise in the force waveforms, though some is still there (about 2 orders of magnitude lower amplitude). I will post plots when the simulations have progressed farther. My guesses as to what this could be due to:
1. Perhaps 12 decimal precision is not enough, though double precision has only 15-17 significant digits of decimal precision.
2. Numerical noise from other calculations.
3. Interface tracking scheme is slightly oscillatory. There are more advanced (and computationally expensive) schemes available.

I did another test where I ran a transient no-motion simulation (interFoam) with the 800k cell mesh until the surface seemed mostly steady. Some motion of the fluid interface, but never more than 1 cell, and no waves/oscillations. See attached plot of the forces. I think those oscillations are probably due to the interface tracking scheme.

Thoughts anybody?

Hello spaceprop
I also face the problem of pressure oscillation in my case (the interaction of dambreak with structure). the phenomenon in my pressure is similar with your picture. But the isoAdvector interface traking (more accurate) is used in my case. And i test the MULES in interfoam,but the interface is not right apparently. I also try differeent schemes in fvscheme, but it does not work.Have you already fixed your problem? or have some idea about this probelm.
DM

[spaceprop]

Hi,

I isolated the remaining oscillations I was seeing (after fixing the precision problem) to either the surface reconstruction scheme or numerical solver. I believe it's fundamentally due to the large density gradient across the free surface. There's no way to get rid of them without changing the solver (maybe higher order and more dissipative, e.g. WENO?) or interface tracking scheme, and even then, I'm not sure it would get rid of them.

That being said, I'm kind of surprised you see them in a 1g dam-break simulation. The liquid motion tends to drown out the small surface oscillations.